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Abstract 

The  purpose  of  this  thesis  was  to  investigate  the  application  of  the  ray 
next-event  estimation  Monte  Carlo  method  to  the  transport  of  primary  and 
secondary  gamma  rays.  The  problem  of  interest  was  estimation  of  the  free  field 
flux  at  a  distant  point  in  a  vacuum  from  a  point  source  in  the  atmosphere.  An 
existing  Fortran  code  for  neutron  transport,  Ray  Next-Event  Estimator  v4.0,  was 
adapted  to  perform  photon  transport  computations  including  coherent  scattering, 
incoherent  scattering,  photoelectric  absorption,  and  pair  production  interactions. 
Production  and  transport  of  secondary  gamma  rays  produced  in  bremsstrahlung, 
neutron  inelastic  scatter,  and  neutron  absorption  interactions  was  also 
implemented. 

A  new  version  of  the  code,  Ray  Next- Event  Estimator  v5.1,  was  produced 
with  the  added  photon  transport  capability  and  other  changes  focused  on  future 
development  of  the  estimator  code  for  application  to  this  class  of  problems.  Code 
version  5.1  was  exercised  and  compared  to  version  4.0  for  neutron  transport 
computations.  Code  version  5.1  was  also  demonstrated  for  application  to  gamma 
ray  transport  computations  and  coupled  neutron-photon  transport  computations. 
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RAY  NEXT-EVENT  ESTIMATOR  TRANSPORT  OF  PRIMARY  AND 

SECONDARY  GAMMA  RAYS 


I.  Introduction 

Monte  Carlo  methods  are  a  popular  choice  to  simulate  the  transport  of 
neutral  particles  in  a  scattering  medium  from  a  source  to  a  point  of  interest.  The 
stochastic  nature  of  neutral  particle  transport  lends  itself  well  to  this  type  of 
simulation  (Zickafoose  &  Mathews,  20I0-20II).  Monte  Carlo  methods  rely  on 
sampling  from  the  domain  of  theoretical  behaviors  to  determine  the  radiation 
field  at  a  specified  location.  To  provide  a  solution  of  sufficiently  low  variance, 
large  numbers  of  particle  histories  must  be  sampled,  resulting  in  correspondingly 
large  computation  times.  Variance  reduction  schemes  can  be  applied  to  reduce 
the  number  of  samples  required  to  converge  on  an  estimate  for  a  particular 
problem. 

The  Air  to  Space  Particle  Transport  Problem 

A  particular  class  of  problems,  estimation  of  the  free-field  flux  at  a  point  in 
a  vacuum,  hereafter  referred  to  as  the  flux-point,  from  a  point  source  in  a 
scattering  medium  is  of  interest.  The  general  problem  is  illustrated  in  Figure  1. 
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Figure  1.  The  Transport  Problem. 


Consider  a  neutral  particle,  a  neutron  or  a  photon,  emitted  from  a  point 
source  in  the  atmosphere.  That  particle  has  energy  and  direction  correlated  with 
the  properties  of  the  source.  The  particle  will  travel  along  its  path  until  it 
interacts  with  atmospheric  nuclei  or  atomic  electrons.  Upon  interaction,  the 
particle  has  a  probability  to  scatter  to  a  flux-point  in  space  with  a  discretized 
grid  in  time  and  energy.  That  probability  is 
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P.  ^  =  probability  of  arrival  at  the  flux  point  in  energy  bin  i  and  time  bin  j 

E'  =  initial  energy 

O'  =  initial  direction 

s  =  distance  traveled  along  the  path  before  collision 

6  =  polar  angle  between  the  particle  path  and  a  vector  to  the  flux-point  at 
the  point  of  scatter 

p  =  the  probability  of  scatter  into  the  angle  0 

T  =  optical  thickness  from  the  point  of  scatter  to  the  flux-point 

E  =  energy  after  scatter 

p  =  distance  from  the  point  of  scatter  to  the  flux-point 
box{^...^  =  a  2-dimensional  normalized  boxcar  function  for  a  time-energy  grid  bin 
with  center  and  width  ^AE,Arj. 

This  is  the  next-event  estimator  (Lux  &  Koblinger,  1991,  p.  380)  and  is  herein 
referred  to  as  the  point  next-event  estimator.  For  a  given  scatter,  there  is  a 
single  probability  and  time-energy  bin  to  which  the  particle  may  contribute.  To 
converge  on  a  solution,  the  estimator  must  be  applied  to  a  significant  number  of 
scattering  points  sufficiently  dispersed  throughout  the  sample  domain  in  energy, 
position,  and  time.  Sampling  of  these  points  is  a  source  of  variance  in  the 
estimate  from  the  point  next-event  estimator. 

Instead,  consider  that  there  is  a  continuous  probability  at  each  point  along 
the  ray,  the  trajectory  of  the  particle,  at  which  the  particle  might  scatter  to  the 
flux-point.  Depending  on  the  scattering  mechanics  of  the  particle  along  the  ray 
there  may  be  a  possibility  of  reaching  multiple  time-energy  bins  at  the  flux-point 
from  a  single  ray  (Zickafoose  &  Mathews,  2010-2011).  Recognizing  that  there  are 
continuous  segments  along  each  ray  from  which  the  particle  could  scatter  to  the 
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flux-point  and  by  treating  each  ray  continuously,  instead  of  only  as  an  endpoint 
as  in  the  point  next-event  estimator,  a  new  estimator  could  eliminate  the  source 
of  variance  from  sampling  the  point  of  scatter. 

The  Ray  Next-Event  Estimator 

The  ray  next-event  estimator  is  a  Monte  Carlo  method  developed  by  Lt  Col 
Matthew  Zickafoose  (2010-2011)  to  improve  the  efficiency  of  existing  tools  and 
methods  for  neutron  transport  computations  from  a  point  source  in  the 
atmosphere  to  a  flux-point  in  a  vacuum.  A  purpose-built  code  accompanied  the 
development  of  the  method  for  demonstration  and  comparison  to  a  conventional 
point  next-event  estimator.  The  estimator  developed  by  Zickafoose  is 

P  - f  ds  ’ -  (2) 

i  (pWf 

where  isotropic  scattering  is  assumed  and 
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P 


probability  of  arrival  at  the  flux-point  in  energy  bin  i  and  time 
bin  i 

initial  energy 
initial  direction 

beginning  of  a  segment  of  the  ray  that  contributes  to  energy 
bin  i  and  time  bin  j 

end  of  a  segment  of  the  ray  that  contributes  to  energy  bin  i 
and  time  bin  j 

optical  thickness  along  the  ray  to  the  beginning  of  the  segment 

macroscopic  scatter  cross  section  at  the  point  s 

optical  thickness  from  5^  to  a  point  s  on  the  segment 

optical  thickness  from  the  point  s  to  the  flux-point 

energy  after  scatter  at  point  s 

distance  from  the  point  s  to  the  flux-point. 


If  isotropic  scattering  cannot  be  assumed,  we  recognize  that  the  factor  of  1  /  dvr 
from  (2)  becomes  /  27r  as  in  the  point  next-event  estimator  from  Lux 

(1991,  p.  380).  The  ray  next-event  estimator  is  then 
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Statement  of  the  Research  Problem 

Incorporate  the  transport  of  primary  and  secondary  gamma  radiations  in 
the  Ray  Next-Event  Estimator  (RNEE)  v4.0  Fortran  code. 


5 


Assumptions  and  Limitations 


The  following  assumptions  and  limitations  were  applied  during  the  course  of 
this  research  to  allow  focus  to  remain  on  the  implementation  of  gamma  ray 
interactions  in  the  ray  next-event  estimator. 

Ray  Next-Event  Estimation 

The  efficacy  of  the  ray  next-event  estimator  was  not  the  subject  of  this 
research.  The  method  was  demonstrated  for  neutron  transport  calculations 
(Zickafoose  &  Mathews,  2010-2011),  but  with  appropriate  application  and 
inclusion  of  additional  physics  data  can  be  generalized  for  neutral  particle 
transport  calculations  including  photon  transport. 

Neutron  Transport  Calculations 

Significant  changes  to  the  implementation  of  neutron  transport  calculations 
were  avoided  throughout  the  research.  The  implementation  of  neutron  transport 
physics  was  assumed  to  be  correct  under  the  assumptions  and  limitations  made 
by  Zickafoose  (2010-2011). 

Atmospheric  Representation 

The  atmosphere  is  assumed  as  the  U.S.  Standard  Atmosphere  and  is 
composed  of  only  nitrogen  and  oxygen  with  relative  concentrations  of  74.8488% 
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and  21.1512%  respectively.  We  assume  a  constant  composition,  variable  density 
atmosphere  up  to  an  altitude  of  86  km.  Altitudes  greater  than  86  km  are 
considered  vacuum. 

Additional  atmospheric  constituents  can  readily  be  included  in  the 
simulation  with  appropriately  processed  physics  data,  but  a  two  element 
approximation  was  sufficient  for  this  research.  Including  atmosphere  at  altitudes 
above  86  km  would  also  be  possible,  with  modification  to  the  code,  by  adding 
additional  atmospheric  layers. 

Geometric  Representation 

The  orientation  of  the  source,  flux-point,  and  earth  are  held  constant  for 
the  duration  of  the  simulation.  This  limits  the  application  of  the  model  to 
problems  in  which  the  flux-point  is  distant  from  the  source,  as  in  the  case  of  a 
source  in  the  atmosphere  and  a  flux-point  in  geosynchronous  orbit,  so  that  their 
relative  motion  is  small  during  the  simulation.  Should  the  code  be  applied  to 
problems  in  which  the  flux- point  is  at  lower  altitudes  than  geosynchronous  orbit, 
this  assumption  may  not  be  valid. 

Photon  Interactions 

Photon  interactions  included  in  the  simulation  are  coherent  scattering, 
incoherent  scattering,  photoelectric  absorption,  and  pair  production.  Other 
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photon  interactions  contribnte  less  than  1%  to  the  total  interaction  cross  section 
for  photons  and  are  neglected  (Lnx  &  Koblinger,  1991,  p.  47). 

For  coherent  scattering,  the  angnlar  distribntion  of  scattered  photons  is 
approximated  as  the  classical  Thompson  cross  section.  The  actual  angular 
distribution  is  the  product  of  the  Thompson  cross  section  and  the  atomic  form 
factor  (Persliden,  1983,  p.  122)  which  corrects  the  angular  distribution  for 
incident  photon  energy  and  the  Z  value  of  the  scattering  material.  The  atomic 
form  factor  becomes  more  important  at  higher  energies  (Persliden,  1983,  p.  122). 
Because  the  coherent  scatter  cross  section  falls  off  rapidly  with  increasing  energy 
and  incoherent  scatter  and  pair  production  interactions  dominate  at  higher 
energies,  the  atomic  form  factor  is  assumed  to  have  a  negligible  impact  on  these 
transport  computations. 

For  incoherent  scattering,  the  Klein-Nishina  approximation  for  the 
distribution  of  scattered  photons  is  assumed  without  the  electron  binding 
correction,  also  called  the  incoherent  scattering  function.  The  incoherent 
scattering  function  corrects  the  distribution  of  scattered  photons  for  differences 
between  elements  and  becomes  important  for  high  values  of  Z  (Persliden,  1983, 
p.  121).  Should  the  model  be  extended  beyond  the  two  element  atmosphere 
approximation,  the  incoherent  scattering  function  may  need  to  be  included. 
Discussion  and  references  for  the  electron  binding  correction  can  be  found  in  Lux 


(1991,  p.  46). 


For  photoelectric  absorption  interactions,  the  emission  of  characteristic 


x-rays  and  anger  electrons  is  ignored.  Shonld  the  estimator  be  applied  to  low 
energy  photon  transport  compntations,  characteristic  x-rays  may  not  be 
negligible.  The  photoelectron  is  considered  in  the  compntation  if  it  is  ejected 
with  snfficient  kinetic  energy  to  emit  bremsstrahlnng  radiation  as  it  decelerates  in 
the  scattering  medinm. 

The  process  of  pair  prodnction  is  approximated  as  the  immediate 
prodnction  of  two  new  gamma  rays  at  the  site  of  the  interaction.  The  path 
length  and  time  of  flight  of  the  electron  and  positron  are  considered  negligible. 

As  in  photoelectric  absorption,  the  electron  and  positron  are  considered  if  they 
are  prodnced  if  they  are  prodnced  with  snfficient  energy  to  emit  bremsstrahlnng 
radiation  as  they  decelerate  in  the  scattering  medinm. 

Secondary  Gamma  Ray  Production  Interactions 

The  model  assumed  from  bremsstrahlnng  radiation  includes  the  continuous 
slowing  down  approximation,  that  all  of  the  kinetic  energy  of  the  decelerating 
particle  is  dissipated,  and  that  the  radiated  energy  is  emitted  as  a  single  coherent 
photon.  These  assumptions  significantly  simplify  the  model  and  limit  its  physical 
relevance,  but  implementation  in  the  code  provides  the  infrastructure  for  a  more 
refined  model  to  be  implemented  in  the  future. 
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The  model  assumed  for  neutron  absorption  assumes  that  the  kinetic  energy 


of  the  absorbed  neutron  is  directly  converted  to  energy  to  place  the  target 
nucleus  in  an  excited  state.  Recoil  energy  of  the  nucleus  is  not  considered.  This 
limits  the  physical  relevance  of  the  model,  but  its  implementation  in  the  code 
provides  the  infrastructure  for  a  more  refined  model  to  be  implemented  in  the 
future. 

Motion  of  Nuclei  and  Charged  Particles 

The  normal  motion  of  atoms  and  the  motion  of  recoiling  nuclei  are 
considered  negligible  compared  to  the  space  and  time  scales  considered  in  the 
problem.  The  path  length  of  charged  particles,  electrons  and  positrons,  is  also 
considered  negligible  when  compared  to  that  of  neutral  particles  in  the  scattering 
medium. 

Decay  of  Excited  Nuclei  by  Gamma  Emission 

Cascades  of  photons  emitted  in  gamma  decay  processes  by  nuclei  excited  in 
neutron  inelastic  scatter  and  other  interactions  are  assumed  to  be  emitted 
instantaneously  after  the  interaction  that  excited  the  nucleus.  This  assumes  that 
the  excited  state  half-lives  and  any  metastable  state  half-lives  are  insignificant 
compared  to  the  time  scale  of  the  simulation. 
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Approach 


The  research  approach  is  tailored  to  allow  the  focus  of  this  research  to 
remain  on  application  of  the  ray  next-event  estimator  to  gamma  ray  transport. 

1.  Adapt  the  purpose  built  Ray  Next-Event  Estimator  4.0  code  for 
gamma  ray  transport. 

2.  Apply  the  Ray  Next-Event  Estimator  code  to  neutron  and  photon 
transport  problems. 

3.  Extend  neutron  and  photon  interaction  physical  models  to  account 
for  generation  of  secondary  gamma  radiation  particles  for  transport 
computations. 

4.  Apply  the  Ray  Next-Event  Estimator  code  to  coupled  neutron- 
photon  transport  problems. 
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II.  Theory 


Four  types  of  interactions  capture  the  vast  majority  of  gamma  ray 
interactions  in  a  scattering  medium  like  air.  In  air,  coherent  scattering  and 
photoelectric  absorption  are  the  dominant  processes  for  photons  with  energy  less 
than  100  keV.  Incoherent  scattering  dominates  photon  interactions  for  energies 
100  keV  to  10  MeV.  Pair  production,  while  not  possible  for  photons  with  energy 
less  than  two  electron  rest  masses  or  1.022  MeV,  quickly  becomes  important 
above  10  MeV.  Cross  sections  for  these  interactions,  as  well  as  the  total  cross 
section,  for  a  nitrogen-oxygen- argon  mixture  that  approximates  air  are  plotted  in 
Figure  2. 


Figure  2.  Photon  Interaction  Cross  Sections  for  a  N-O-Ar  Mixture  Like  Air  (NIST,  2010). 
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Other  particle  interactions  can  also  prodnce  gamma  ray  photons  as  a  result 


of  the  interaction.  High  energy  electrons  and  positrons  produced  in  other  gamma 
interactions  can  produce  bremsstrahlung  radiation  photons  as  they  decelerate  in 
the  scattering  medium.  Gamma  rays  can  also  be  emitted  by  nuclei  left  in  an 
excited  state  from  neutron  inelastic  scatter  and  neutron  absorption  interactions. 

Coherent  Scattering 

Coherent  scattering,  also  called  Rayleigh  scattering,  occurs  when  a  photon 
interacts  coherently  with  the  electrons  of  a  target  atom.  The  photon  is  absorbed 
and  immediately  reemitted  from  the  atom  in  a  different  direction  but  with  the 
same  energy  as  the  incident  photon  (Choppin,  Liljenzin,  &  Rydberg,  2002,  p. 

143).  Figure  3  shows  a  graphical  depiction  of  a  coherent  scatter  interaction. 
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The  angular  distribution  is  described  as  the  product  of  the  classical 


Thompson  cross  section  and  the  atomic  form  factor, 

^  =  ±(l  +  co.^g)F(q,z)  ,  (4) 

where  r  is  the  classical  electron  radius,  and  zj  is  the  atomic  form  factor 
(Lux  &  Koblinger,  1991,  p.  47). 

The  atomic  form  factor  is  dependent  on  the  photon  energy  and  the  atomic 
number  of  the  target  atom.  Including  the  atomic  form  factor  in  these 
computations  would  significantly  complicate  the  model  for  coherent  scattering. 
We  assume  the  effect  of  the  atomic  form  factor  to  be  negligible  (T^g, zj  =  1)  to 

allow  the  focus  of  the  research  to  remain  on  implementation  of  gamma  ray 
interactions  in  the  code.  The  atomic  form  factor  becomes  more  important  for 
higher  energy  photons  (Persliden,  1983,  p.  122),  but  because  the  cross  section  for 
coherent  scatter  drops  rapidly  with  increasing  energy  relative  to  other  photon 
interactions  the  probability  of  a  high  energy  coherent  scatter  interaction  is  low. 
References  for  formulations  of  the  atomic  form  factor  can  be  found  in  Lux  (1991, 
p.  46). 

The  angular  distribution  of  coherently  scattered  photons,  scaled  so  that 
0  <  <  1 ,  is  plotted  in  Figure  4  for  photons  incident  from  the  left. 
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This  is  the  angular  cross  section  for  a  differential  element  of  polar  scattering 
angle. 

Computationally,  it  is  more  appropriate  to  consider  (5)  in  terms  of  the 
cosine  of  the  scattering  angle,  because  all  of  the  subsequent  transport  calculations 
use  the  cosine  and  not  the  angle  itself.  Rearranging  (5)  we  get 

d(j  =  Trr  |l  +  cos^  djsinddd  .  (6) 


If  we  let  n  =  cos  9  ,  then  dp. 


=  sin9d9  and  dO  = 


sind 


dp  ,  and  (6)  becomes 


da  =  Trr  |l  +  p^  jdp 


The  angular  cross  section  for  a  differential  element  of  p  is  then 


da  r 


=  ^(l  +  p'] 
dp  TT  '  ^ 


(7) 


(8) 


Scaling  so  that  0  <  —  <  1  and  defining  the  scaled  distribution  function, 
dp  dp 


/(p),  gives 


(9) 


The  distribution  /(p)  is  plotted  in  Figure  5.  The  scaling  step  from  (8)  to  (9)  is 
not  entirely  necessary  for  the  implementation  of  coherent  scatter  angular 
selection.  However,  scaling  the  angular  distribution  function  will  be  convenient 
for  the  implementation  of  incoherent  scatter  scattering  angle  selection  later. 
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Forgoing  the  change  from  (8)  to  (9)  and  defining  =  —  =  -  — +  does 

^  '  dfX  TT  ''  ' 


not  change  the  resnlts  of  the  implementation  for  coherent  scatter  scattering  angle 


selection. 


/(y^) 


Figure  5.  Distribution  of  the  Cosine  of  the  Scattering  Angle  for  Coherently  Scattered 

Photons. 


Photoelectric  Absorption 

Photoelectric  absorption  occnrs  when  a  gamma  ray  photon  is  fnlly  absorbed 
by  an  absorber  atom  and  canses  the  ejection  of  an  electron,  called  a 
photoelectron,  from  a  bonnd  shell  of  the  atom.  The  kinetic  energy,  T  ,  of  the 
photoelectron  is 
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T  =  hv  -  E, 

0 


(10) 


where  is  the  binding  energy  of  the  photoelectron  in  its  original  bonnd  shell 
(Knoll,  2000,  p.  309).  The  resnlting  ionized  absorber  atom  may  also  emit  one  or 
more  characteristic  x-ray  photons  or  an  Auger  electron  as  electrons  from  other 
bound  shells  rearrange  to  fill  the  vacancy  left  by  the  ejected  photoelectron 
(Knoll,  2000,  p.  309).  Figure  6  depicts  the  photoelectric  absorption  process. 


Figure  6.  Photoelectric  Absorption. 


For  the  purposes  of  this  research,  we  assume  that  the  angular  distribution 
of  photoelectrons  is  isotropic  and  that  characteristic  x-rays  and  auger  electrons 
are  negligible.  These  assumptions  simplify  the  model  for  photoelectric  absorption 
to  allow  the  focus  of  the  research  or  remain  on  implementation  of  gamma  ray 
interactions  in  the  code. 
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Incoherent  Scattering 


Incoherent  scattering,  also  known  as  Compton  scattering,  occnrs  when  an 
incident  photon  interacts  with  an  electron  in  the  scattering  medinm.  The 
incident  photon  transfers  some  of  its  energy  to  the  electron  and  is  scattered 
throngh  an  angle  9  (Knoll,  2000,  p.  310)  as  shown  in  Fignre  7. 


hv 

WWWWXAy^ 


Before 


After 


Figure  7.  Incoherent  Scattering. 


The  energy  of  the  scattered  photon  is 


hv'  = 


hv 


1  +  a  ^1  -  cosdj 


(11) 


where  hv'  is  the  energy  of  the  scattered  photon,  hv  is  the  energy  of  the  incident 
photon,  and  a  =  hv  j  ;  the  ratio  of  incident  photon  energy  to  the  rest  mass 
energy  of  the  electron  (Knoll,  2000,  p.  310).  The  kinetic  energy,  T,  of  the  recoil 
electron  is 


T  =  hv  -  hv' 


(12) 
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The  angular  distribution  of  scattered  photons  is  approximated  by  the 


Klein-Nishina  formula, 


dn 


^  1  +  cos^  6  ' 


1  +  a  (1  - 


^1  -  cosdj 


1  + 


“M 

1  -  cosdj 

|l  +  cos^  dj 

1  +  a^l  -  cosdj 

,(13) 


where  Z  is  the  atomic  number  of  the  scattering  medium  (Knoll,  2000,  p.  51). 
The  angular  distribution  of  incoherently  scattered  photons,  scaled  so  that 

0  <  <  1 ,  is  plotted  in  Figure  8  for  photons  incident  from  the  left  and  several 


values  of  a .  Recalling  that  a  =  hv  /  ,  so  that  larger  values  for  a 

correspond  to  higher  incident  photon  energies,  note  the  strong  tendency  for 
forward  scattering  for  higher  photon  energies. 

This  model  for  incoherent  scattering  assumes  that  the  electron  with  which 
the  incident  photon  interacts  is  initially  unbound  and  at  rest.  In  reality,  most 
electrons  in  a  scattering  medium  meet  neither  of  these  criteria.  It  is  possible  to 
refine  the  Klein-Nishina  approximation  using  electron  binding  corrections  or 
incoherent  scattering  functions,  but  generally  unnecessary  for  practical  transport 
computations  (Lux  &  Koblinger,  1991,  p.  46).  Electron  binding  corrections  and 
incoherent  scattering  functions  are  not  considered  here  in  order  to  allow  the  focus 
of  the  research  to  remain  on  the  implementation  of  gamma  ray  transport  in  the 
code. 
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^  =  2^Zr^ 

de 


1  +  a{l-  cos 9^ 
f 

1  +  -, - 


1  +  cos^  9 


1  +  cos  9 


a~  ^1  -  cos^j 

y  1  +  a  ^1 


sin0  . 


1  -  cos  ( 


This  is  the  angular  cross  section  for  a  differential  element  of  polar  scattering 
angle. 

As  before,  it  is  computationally  more  appropriate  to  consider  (14)  in  terms 
of  the  cosine  of  the  scattering  angle,  because  all  of  the  subsequent  transport 
calculations  use  the  cosine  and  not  the  angle  itself.  Rearranging  (14)  gives 


da  =  27iZr 


1  +  a^l  -  cos^j 


1  +  cos^  9 


(l  -  cos^) 

1  +  ^  sin  9d9  . 

1^1  +  cos^  0)  l  +  Q;(l-cos0j 


If  we  let  p.  =  cos  9  ,  then  \dw  =  sm9d9  and  d9  = - dp  ,  and  (15)  becomes 

'  '  sin  9 '  ' 


Zr^  1 
da  =  — ^  - - - r 

27r  I  1  +  «(!  -  pj 


1  +  p 
2 


o? 

(i-p; 

) 

|l  _l_  p2  j 

1  +  cr  1 

1-p) 

r 

dp  .  (16) 


Angular  cross  section  for  a  differential  element  of  p  is  then 


da  Zr: 


dp  27r  M  +  «(l  -  pj 


1  +  p 
2 


o? 

(l-p 

) 

(l  +  p) 

1  +  1 

i-p) 
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Scaling 


da 

dfi 


so  that  0  < 


—  <  1  and  defining  the  scaled  distribntion  fnnction, 
d^ 


gi 


gives 


1  +  a{l- 


yv  2  , 


1  + 


«■  ^1  -  /wj 


{l  +  1  +  a  ^1  - 


(18) 


The  distribntion  is  plotted  for  several  valnes  of  a  in  Fignre  9. 


f  (fi  ■,  a) 


Figure  9.  Distribution  of  the  Cosine  of  the  Scattering  Angle  for  Incoherently  Scattered 

Photons  of  Various  Energies. 
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Pair  Production 


In  pair  production,  a  gamma  ray  interacts  within  the  coulomb  field  of  a 
nucleus  or  the  atomic  electrons  and  transitions  into  an  electron-positron  pair 
(Knoll,  2000,  p.  312).  To  create  an  electron-positron  pair,  the  incident  photon 
must  have  energy  greater  than  two  electron  rest  masses  or  1.022  MeV.  Any 
energy  in  excess  of  1.022  MeV  is  converted  to  kinetic  energy  shared  by  the 
resulting  electron-positron  pair  (Knoll,  2000,  p.  312). 

The  positron  created  in  the  pair  production  interaction  subsequently 
decelerates  and  annihilates  with  an  electron  in  the  scattering  medium  releasing 
the  mass  energy  of  the  positron  and  electron  in  the  form  of  two  511  keV  gamma 
photons  (Knoll,  2000,  p.  52).  These  secondary  photons  have  opposite  directions 
of  travel  that  are  uncorrelated  to  the  direction  of  the  original  photon  (Lux  & 
Koblinger,  1991,  p.  47).  Figure  10  depicts  a  pair  production  interaction  and 
subsequent  annihilation. 
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Pair  production  is  in  fact  a  much  more  complicated  process,  but  this  model 
is  adequate  for  practical  transport  calculations.  A  more  detailed  description  can 
be  found  in  the  article  by  Hubbell  (1980). 


Bremsstrahlung  Radiation 


Photoelectric  absorption,  incoherent  scattering,  and  pair  production 
interactions  have  potential  to  produce  high  energy  electrons.  As  an  electron 
moves  through  a  scattering  medium,  it  can  be  deflected  through  large  angles  by 
collisions  with  other  electrons.  Because  charged  particles  radiate  electromagnetic 
energy  when  accelerating,  the  decelerating  electron  emits  radiation  along  its  path 
(Krane,  1988,  p.  196).  This  radiation  is  called  bremsstrahlung,  German  for 
braking,  radiation.  A  simplified  model  for  production  of  photons  by 
bremsstrahlung  interactions  is  presented  here  and  implemented  in  the  code.  The 
model  assumes  the  continuous  slowing  down  approximation  for  the  charged 
particle  and  that  the  radiated  energy  is  emitted  as  a  single  coherent  photon. 

The  energy  radiated  per  unit  path  length  is  given  by 


(19) 


where  is  the  fine  structure  constant,  N  is  the  number  density  of  the  scattering 
medium,  and  E  is  the  energy  of  the  incident  electron  (Knoll,  2000,  p.  44). 
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Electrons  also  dissipate  energy  along  their  path  as  they  collide  with  other 
particles  in  the  scattering  medinm.  Energy  loss  per  nnit  path  length  dne  to 
collisions  is 


where  {3  =  v  /  c  and  /  is  the  empirically  determined  electron  excitation  and 
ionization  potential  of  the  target  atom  (Knoll,  2000,  p.  44). 

The  total  energy  loss  per  nnit  path  length,  or  stopping  power  of  the 
scattering  medinm,  is 

dx  y  dx  dx 

Determining  the  range  of  an  electron  in  a  scattering  medinm  can  be 
accomplished  by  integrating  both  the  radiative  and  collisional  energy  losses  over 
the  path  of  the  electron.  As  a  result  of  the  erratic  nature  of  the  path  of  the 
scattered  electron,  this  is  a  difficult  process  (Krane,  1988,  p.  197).  Instead,  range 
to  incident  energy  relationships  can  be  estimated  from  empirical  data  (Krane, 
1988,  pp.  197-198). 

Bremsstrahlung  radiation  is  small  for  electrons  traveling  at  non-relativistic 
speeds  or  energy  below  1  MeV  (Krane,  1988,  p.  196).  For  relativistic  electrons, 
the  ratio  of  radiative  losses  to  collisional  losses  can  be  estimated  by 
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(22) 


fdE^ 

{dx)^ 

where  T  is  the  kinetic  energy  of  the  incident  electron  (Krane,  1988,  p.  197). 

Decay  of  Excited  Nuclei  by  Gamma  Emission 

The  atomic  nnclens  can  be  excited  to  higher  energy  states,  or  levels,  in 
neutron  and  other  interactions.  These  excited  states  decay  by  a  number  of 
processes  including  the  emission  of  gamma  rays.  A  simplified  model  of  the 
gamma  emission  process  is  presented  here  that  is  sufficient  for  this  research,  but 
is  in  no  way  a  complete  model  for  the  process.  A  more  complete  discussion  of 
gamma  decay  of  nuclear  excited  states  may  be  found  in  Krane  (1988). 

Gamma  decay  modes  of  nuclear  excited  states  are  different  for  each  element 
and  each  isotope  of  that  element.  A  graphical  depiction  of  the  different  energy 
levels  is  often  referred  to  as  a  level  scheme.  The  level  scheme  for  is  shown  in 
Figure  11.  Each  excited  state,  or  level,  is  labeled  with  its  energy  in  keV,  and 
downward  arrows  indicate  gamma  emissions  as  the  excited  nucleus  transitions 
from  one  excited  state  to  another. 

Properties  of  the  nucleus  and  quantum  mechanical  principles  prohibit 
particular  energy  transitions  and  determine  half-lives  between  transitions,  and 
there  are  many  models  to  explain  emission  behavior.  In  short,  there  are  multiple 


dx 


T  +  Z 
1600 
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emission  chains  that  can  be  followed  from  one  excited  state  down  to  the  lowest 


energy  state  available.  Each  of  these  chains  has  a  relative  intensity  or  likelihood 
of  occurring.  Empirical  data  on  the  decay  of  nuclear  excited  states  is  available  in 
tabular  format  from  the  Nudat  2.5  database  (NNDC,  2010). 
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Figure  11.  Level  Scheme  and  Decay  Gammas  for  energy  levels  in  keV.  (NNDC,  2010) 
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III.  Methodology 


The  Ray  Next-Event  Estimator  (RNEE)  v4.0.3  code  was  written  in  Fortran 
2003  nsing  the  Microsoft  Visnal  Stndio  2005  development  environment  and  the 
Intel  Fortran  compiler  v9. 1.023.  The  final  version  prodnced  in  the  conrse  of  this 
research  was  RNEE  v5.1.05,  developed  and  compiled  in  Microsoft  Visnal  Stndio 
2008  nsing  the  Intel  Fortran  compiler  vl2.0.2.154. 

RNEE  Transport  of  Speed-of-Light  Neutral  Particles 

The  existing  method  for  ray  estimation  is  well  snited  for  adaptation  to  the 
pnrpose  of  gamma  ray  transport  calcnlations.  Compntational  transport  of  speed 
of  light  nentral  particles  is  largely  similar  to  compntational  nentron  transport. 
The  two  fnndamental  differences  between  nentron  and  photon  transport  are  that 
photon  velocity  is  independent  of  energy  and  always  eqnal  to  the  speed  of  light, 
and  that  photon  interactions  with  matter  follow  different  physical  models  than 
nentron  interactions. 

Implementation  of  Gamma  Ray  Interaction  Physics 

The  physical  models  of  gamma  interactions  are  nsed  by  the  code  for  three 
main  compntations: 
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1.  Drawing  the  scattering  angle  from  the  appropriate  angnlar 
distribntion  for  instances  where  the  scattering  angle  is  nnknown. 

2.  Determination  of  the  likelihood  of  scattering,  angnlar  density,  into  a 
particular  scattering  angle  when  that  angle  is  known. 

3.  Determination  of  the  particle  energy  after  a  scattering  interaction 
when  the  scattering  angle  is  known. 

Scattering  Angle  Selection 

The  dependence  of  the  scattering  angle  distributions  on  the  energy  of  the 
incident  particle  is  one  of  the  primary  differences  between  the  transport  of 
photons  and  isotropically  scattered  neutrons  as  implemented  in  RNEE  v4.0.  In 
the  three-dimensional  representation  of  the  photon  interaction,  the  scattering 
angle  is  defined  by  two  angles;  9  ,  the  polar  scattering  angle,  and  a; ,  the 
azimuthal  scattering  angle.  The  polar  scattering  angle  is  the  angle  between  the 
incident  and  departure  photon  velocity  vectors  in  the  plane  defined  by  those 
vectors.  The  azimuthal  scattering  angle  is  the  angle  of  rotation  about  the 
incident  photon  velocity  vector.  In  general,  6  is  the  angle  to  which  the  energy 
dependent  angular  distributions  apply  and  is  drawn  from  these  distributions  on 
the  interval  zero  to  tt  .  The  azimuthal  angle,  co  ,  is  drawn  uniformly  from  the 
interval  zero  to  27r .  Although  properties  of  the  gamma  ray  and  the  scattering 
medium,  such  as  polarization,  may  impact  the  distribution  of  co  ,  uniform 
selection  of  this  angle  is  appropriate  for  these  general  transport  calculations. 
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Coherent  Scattering 


Recall  the  scaled  distribution  function  for  coherently  scattered  photons  in 
terms  of  the  cosine  of  the  scattering  angle: 

,  (i  +  m" 

We  convert  /(p-)  to  a  probability  density  function  (PDF),  ,  by 


P 


(23) 


The  PDF,  ,  can  be  converted  to  a  cumulative  probability  density  function 
(CDF),  P(A^),by 


P(^)  = 


(24) 


where  fi  is  a  dummy  variable  of  integration.  This  gives 

P(p.)  =  — ^4  +  3p,  +  j  ■  (25) 

Figure  12  shows  plots  of  the  PDF  and  CDF  for  the  distribution  of  the  cosine  of 
the  scattering  angle  for  coherent  scattering. 
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Figure  12.  PDF  and  CDF  for  the  Distribution  of  the  Cosine  of  the  Polar  Scattering  Angle  for 

Coherently  Scattered  Photons. 

To  sample  p. ,  we  draw  a  random  number,  ^  ,  uniformly  distributed  between 
zero  and  one  and  evaluate  to  find  p. .  Evaluation  of  can  be 

accomplished  by  root-solving  or,  if  P(/wj  is  invertible,  direct  computation.  The 
P(^p.j  in  (25)  can  be  directly  inverted  and  has  three  roots,  one  real  and  two 

imaginary.  Because  the  imaginary  roots  represent  non-physical  values  for  /x  , 
only  the  real  root  is  relevant: 

l-(2-4C  +  V5  +  16^K-l)P 

^  =  P  (d  =  V - ,  1  ■  (26) 

(2-4{  +  ,/6  +  16eK-l)) 

This  allows  us  to  compute  appropriately  distributed  cosines  for  the  scattering 
angle  of  coherently  scattered  photons.  Figure  13  shows  the  cumulative 
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distribution  of  100,000  angle  cosines  drawn  by  the  estimator  code  along  with  a 


line  representing  the  cumulative  probability  density  function. 


P{^i) 


Figure  13.  Distribution  of  100,000  Angle  Cosines  Drawn  by  the  Estimator  Code  Compared  to 
the  Angle  Cosine  CDF  for  Coherent  Scattering. 


To  verify  that  the  angle  cosines  were  drawn  from  the  correct  distribution, 
we  apply  a  statistical  test  for  goodness  of  fit.  We  choose  the  Chi-squared  ( x  ) 
test  for  simplicity  and  readily  available  tables  of  critical  values.  The  x  test 
statistic  is 


X 


2 


(27) 
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where  k  is  the  number  of  bins  into  which  the  data  are  sorted,  O  is  the  number 
of  samples  in  the  bin  i ,  N  is  the  total  number  of  samples,  and  F.  is  the  expected 
fraction  of  samples  in  bin  i .  The  null  hypothesis  is  that  the  data  are  from  the 
distribution  F  .  The  alternate  hypothesis  is  that  the  data  are  from  a  distribution 
other  than  F  .  If  the  test  statistic  is  greater  than  the  upper  critical  value  or  less 
than  the  lower  critical  value  for  the  desired  level  of  confidence  and  number  of 
degrees  of  freedom,  we  reject  the  null  hypothesis  with  the  specified  level  of 
confidence. 

The  100,000  angle  cosines  drawn  by  the  code  were  sorted  into  50  equally 
spaced  bins,  the  f^st  was  applied,  and  we  accepted  the  null  hypothesis.  The 
results  of  the  test  are  shown  in  Table  1. 


Table  1.  Results  of  Chi-Squared  Test  for  Fit  of  100,000  Angle  Cosines  Drawn  by  the  Code 

for  Coherent  Scattering. 


x' 

Critical  Values  (90%  confidence,  49  DoF) 

Accept  Null 

Hypothesis? 

Upper 

Lower 

48.64 

66.34 

33.93 

Yes 

Incoherent  Scattering 


Recall  the  scaled  distribution  function  for  incoherently  scattered  photons  in 
terms  of  the  cosine  of  the  scattering  angle: 
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N  1  1  +  ^1  a 

=  - 7 - r  -  1  +  ^ - TF - ^ - FT  ■ 

tl  +  «(l-^)Jl  ^  (l  +  ^^)[l  +  «(l-A^)]^ 

Converting  to  a  PDF  and  CDF  is  complicated  and  the  resnlting 

fnnctions  are  not  easily  conditioned  for  compntational  applications.  The 
resnlting  CDF  mnst  be  recomp nted  for  a  new  valne  of  a  for  each  interaction  and 
is  also  not  easily  inverted,  increasing  the  compntational  cost  and  necessitating  an 
iterative  root-solving  method  to  sample  p. . 

To  avoid  complicated  implementation  and  the  high  compnting  cost 
associated  with  attempting  to  evalnate  for  the  Klein-Nishina  distribntion, 

a  simpler  approach  is  to  nse  a  method  called  sampling  by  rejection.  Sampling  by 
rejection  trades  the  cost  of  evalnating  P~^  for  the  cost  of  generating  extra 
random  nnmbers  (Lewis  &  Miller,  1993,  p.  303). 

To  sample  from  we  draw  a  pair  of  random  nnmbers  each  nniformly 

distribnted  from  zero  to  one,  and  .  We  then  scale  and  to  represent 
a  point,  ,  in  the  sample  space  for  ;  -1  <  <  1  and  0  <  ^2  -  1  •  We 

then  evalnate  at  ,  and  if  ^2  -  then  we  accept  as  a  valne  for 

l^-  If  ^2  >  f(Ci  ;Q;j ,  we  reject  and  draw  a  new  random  point  from  the  sample 
space  to  test.  Fignre  14  illnstrates  the  sampling  by  rejection  process  for  two 
points  from  /(p;Q;  =  0.2j .  Points  drawn  in  the  shaded  region  are  rejected. 
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f(ix;a^  0.2) 


Figure  14.  Sampling  by  Rejection  from  the  Klein-Nishina  Distribution. 


The  sampling  efficiency,  ry ,  is  the  ratio  of  the  area  of  the  sample  space  to 
the  area  below  the  distribution  curve.  Recalling  from  Figure  9  that  the  shape  of 
the  distribution  changes  with  a ,  we  can  define  the  sampling  efficiency 


Figure  15  shows  a  plot  of  sampling  efficiency  as  a  function  of  a . 


(28) 
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Figure  15.  Sampling  Efficiency  for  Basic  Sampling  by  Rejection  from  the  Klein-Nishina 

Distribution. 


Graphical  inspection  of  immediately  snggests  that  the  maximnm  efficiency 

is  significantly  less  than  nnity.  In  fact,  lim  rylaj  =  2  /  3:  The  maximnm 

efficiency  for  sampling  by  rejection  nsing  this  method  is  66%  for  only  the  lowest 
energy  photons.  As  energy  of  the  incident  photon  increases,  the  sampling 
efficiency  falls  rapidly. 

A  simple  improvement  to  the  rejection  sampling  method  ontlined  above  is 
to  exclnde  from  the  sample  space  any  points  that  are  guaranteed  to  be  rejected. 
Consider  the  limit  of  as  a  approaches  zero: 
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L  2\ 

g{fi)  =  lim/(^;a)  = - - - 

'  >  a^O  '  ^  2 


(29) 


Any  points  selected  in  the  sample  space  above  this  curve  will  be  rejected  for  any 
distribution  with  positive  a .  We  can  eliminate  these  rejections  by  bounding  the 
sample  space  to  exclude  these  samples  from  the  beginning  of  the  sampling 
process. 

Conveniently,  (29)  is  the  same  distribution  as  previously  discussed  for  the 
coherent  scatter  angular  distribution,  (9).  Drawing  and  each  uniformly 
between  zero  and  one,  we  draw  from  the  distribution  by  evaluating  (26) 

at  .  The  second  random  number,  ^2^  is  then  found  by  multiplying  by 
5'(^i)  fo  scale  the  sample  height  below  the  sample  space  bounding  function. 

Figure  16  illustrates  sampling  by  rejection  in  the  new  sample  space.  Note 
the  significantly  smaller  rejection  region  (shaded)  than  in  Figure  14. 
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Figure  16.  Sampling  by  Rejection  from  Bounded  Sample  Space  from  the  Klein-Nishina 

Distribution. 


Sampling  efficiency  for  the  method  using  the  bounded  sample  space  is 


(30) 


The  sampling  efficiency  for  sampling  by  rejection  from  the  bounded  sample  space 
is  plotted  in  Figure  17. 
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Figure  17.  Sampling  Efficiency  for  Improved  Sampling  by  Rejection  from  the  Klein-Nishina 

Distribution. 


While  the  maximum  efficiency,  limryla) ,  is  now  equal  to  unity,  this  method  still 
quickly  becomes  inefficient  for  even  moderately  large  values  of  a . 

The  function  used  to  bound  the  sample  space,  ,  need  not  be 
independent  of  a .  In  fact,  a  boundary  function  that  would  adjust  its  shape 
according  to  a  could  greatly  improve  sampling  efficiency  for  higher  values  of  a . 
The  boundary  function  must  meet  two  requirements:  ^  /  (h-)  ^r  all  p.  in 

the  sampling  interval,  and  it  must  be  possible  to  sample  from  the  distribution 
Preferably,  giyl^  has  a  form  of  its  CDF  that  is  invertible  and  can  be 
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evaluated  with  little  computational  cost,  so  that  the  cost  of  sampling  from  gyl^j 
will  not  diminish  the  improvement  in  efficiency  to  the  sampling  by  rejection 
process.  The  boundary  function  =  lim/^«,p,j  provides  sufficient 

improvement  to  the  sampling  by  rejection  process  for  this  research. 

Table  2  compares  the  observed  rejection  sampling  efficiency  for  the  100,000 
selected  angle  cosines  compared  to  the  theoretical  sampling  efficiency  from  (30) 
and  Figure  17.  Figure  18  shows  the  cumulative  distribution  100,000  incoherent 
scattering  angle  cosines  drawn  by  the  estimator  code  for  a  =  0.001,  0.01,  0.1,  1, 
10,  and  100  respectively  with  lines  representing  the  cumulative  probability 
density  functions. 


Table  2.  Comparison  of  Theoretical  and  Observed  Rejection  Sampling  Efficiencies  for  100,000 
Angle  Cosines  Drawn  by  the  Estimator  Code  for  Various  Values  of  Alpha  for  Incoherent 

Scattering. 


a 

rj  -  Theoretical 

rj  -  Observed 

OL 

rj  -  Theoretical 

rj  -  Observed 

0.001 

0.9980 

0.9977 

1 

0.4307 

0.4306 

0.01 

0.9805 

0.9811 

10 

0.1228 

0.1225 

0.1 

0.8413 

0.8409 

100 

0.0215 

0.0214 
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P  (jx  ]  a  ^  omi)  P  ifx  ■,  a  ^  0.01) 
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Figure  18.  Distribution  of  100,000  Angle  Cosines  Drawn  by  Estimator  Code  for  Various 
Values  of  Alpha  Compared  to  the  Angle  Cosine  CDFs  for  Incoherent  Scattering. 


The  test  was  also  applied  to  each  set  of  angle  cosines  drawn  by  the  code 
to  test  the  goodness  of  fit  to  .  For  each  value  of  a  we  accepted  the  null 

hypothesis.  Table  3  shows  the  results  of  the  tests. 
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Table  3.  Results  of  Chi-Squared  Test  for  Fit  of  100,000  Incoherent  Scatter  Angle  Cosines 
Drawn  by  the  Code  for  Various  Values  of  Alpha. 


OL 

x' 

Critical  Values  (90%  confidence,  49  DoF) 

Accept  Null 

Hypothesis? 

Upper 

Lower 

0.001 

42.47 

66.34 

33.93 

Yes 

0.01 

52.08 

Yes 

0.1 

41.27 

Yes 

1 

46.39 

Yes 

10 

39.90 

Yes 

100 

60.23 

Yes 

Isotropic  Angles 


Angle  cosines  for  other  photon  interaction  processes  are  drawn  uniformly 
from  negative  one  to  one.  This  includes  nuclear  excited  state  gamma  emission, 
emission  of  bremsstrahlung  photons,  and  the  first  annihilation  photon  emitted 
after  a  pair  production  interaction.  The  second  annihilation  photon  is 
constrained  to  travel  in  the  opposite  direction  from  the  first.  Thus,  for  the 
second  photon 


and 


^hi>2 


(31) 


U 


hl’n 


=  UJ 


hi/. 


-I-  TT 


(32) 
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This  selects  the  direction  of  travel  for  the  second  photon  as  the  direction  of  travel 


of  the  first  photon  reflected  throngh  the  origin,  aligning  the  two  annihilation 
photons  correctly  for  the  pair  prodnction  model. 


Angular  Density 


Recall  the  ray  next-event  estimator  from  (3) 


P 

h] 


(A\E',n',s) 


2ti 


We  recognize  that  for  anisotropically  scattered  particles,  the  probability  of 
scatter  into  a  particular  angle,  must  be  computed  to  evaluate  the  integral. 

This  probability,  hereafter  referred  to  as  angular  density,  requires  a  different 
computational  approach  for  each  type  of  photon  interaction. 


Coherent  Scattering 

To  compute  the  angular  density  for  a  coherent  scatter  interaction,  we  need 
only  evaluate  the  probability  density  function  for  the  angular  distribution  of 
photons  from  coherent  scattering  for  the  known  value  of  /i .  The  coherent  scatter 
PDF,  (23),  is 

PiA)  =  + 
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Evaluating  this  function  for  the  known  value  of  /x  for  scatter  to  the  flux-point 
will  give  the  angular  density  for  coherent  scatter  for  that  particular  value  of  p. . 


Incoherent  Scattering 


Computing  the  angular  density  of  an  incoherent  scatter  contribution 
encounters  the  same  challenges  as  drawing  an  angle  cosine  from  the 
Klein-Nishina  distribution.  Recall  the  scaled  version  of  the  Klein-Nishina 


distribution,  (18): 


1  +  a{l  - 


1  +  p. 

V  ^  y 


1  + 


a  (l  “  h-) 


|l  +  j  1  +  a{l  - 


We  can  convert  to  a  PDF  by 


P 


(^;a) 


(33) 


However,  the  algebraic  form  of  this  function  is  complex  and  poorly  conditioned 
for  numerical  implementation.  Instead,  we  evaluate  the  numerator  and 
denominator  in  the  right  hand  side  of  (33)  independently  and  then  divide. 

The  numerator  is  easily  evaluated  for  the  known  values  of  a  and  fx .  The 
denominator,  carrying  out  the  integration,  is 


f‘  +  +  ,  ln(2a  +  l)(«(«-2)-^) 

J-CW’  2  ^  ^^3 

y2a  +  Ij 


(34) 
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This  is  easily  evaluated  for  most  values  of  a ,  but  becomes  numerically  unstable 


for  small  a.  Therefore,  for  a  greater  than  0.01,  we  evaluate  the  right  hand  side 
of  (34).  For  a  less  than  0.01,  we  implement  a  numerical  quadrature  to  evaluate 
the  left  hand  side  of  (34).  Once  values  are  determined  for  the  numerator  and 
denominator  on  the  right  hand  side  of  (33),  is  found  by  carrying  out  the 

division.  This  gives  the  angular  density  for  incoherent  scatter  for  that  particular 
value  of  p. . 

Pair  Production 

The  annihilation  photons  produced  in  pair  production  interactions  also 
require  an  angular  density  adjustment,  but  for  a  different  reason  from  the 
scattering  interactions.  The  computed  contribution  to  a  time-energy  bin  at  the 
flux- point  is  normally  based  on  the  likelihood  of  a  scatter  to  the  flux- point.  In 
the  case  of  pair  production,  the  first  annihilation  photon  is  treated  as  an 
isotropically  scattered  photon  and  the  contribution  to  the  flux-point  is  computed. 
However,  we  must  also  consider  the  case  in  which  the  first  annihilation  photon 
travels  away  from  the  flux- point,  aligning  the  second  annihilation  photon  to 
contribute.  Instead  of  a  single  angle  at  which  the  interaction  can  contribute  to 
the  flux-point,  there  are  two.  This  reduces  the  solid  angle  over  which  the  scatter 
is  distributed  from  dyr  to  27r .  Assuming  that  the  two  photons  produced  in  the 
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pair  production  interaction  are  identical  except  for  their  opposite  directions  of 
travel,  =  1 . 

Post- Scatter  Energy 

The  estimator  code  must  also  compute  the  energy  of  a  photon  after  a 
simulated  interaction.  For  these  computations,  the  energy  of  the  incident  photon 
is  known,  as  well  as  the  cosine  of  the  departure  angle  from  the  interaction. 


Incoherent  Scattering 


Photon  energy  after  an  incoherent  scatter  interaction  is  computed  by 
evaluating  (11)  for  hu'  using  the  energy  of  the  incident  photon,  hu  ,  and  the 
cosine  of  the  known  scattering  angle: 


hu'  = 


Constant  Energy  Processes 

For  coherent  scattering  and  pair  production,  the  post  scatter  energy  is  not 
dependent  on  the  scattering  angle.  The  energy  of  a  coherently  scattered  photon 
is  identical  to  the  energy  of  that  photon  before  the  scattering  interaction 
(Choppin,  Liljenzin,  &  Rydberg,  2002,  p.  143).  For  pair  production,  the  energy 
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for  each  of  the  annihilation  photons  is  eqnal  to  the  rest  mass  of  an  electron  or 
511  keV  (Knoll,  2000,  p.  52). 

Secondary  Particles 

Secondary  or  induced  particles  are  particles  created  as  a  result  of  an 
interaction  simulated  by  the  estimator  code. 

Bremsstrahlung  Photons 

Pair  production,  photoelectric  absorption,  and  incoherent  scattering 
interactions  may  produce  high  energy  charged  particles,  electrons  and  positrons, 
which  will  decelerate  in  the  scattering  medium.  The  energy  radiated  by  a 
decelerating  charged  particle,  bremsstrahlung  radiation,  can  be  computed  by 
combining  (21)  and  (22)  and  assuming  that  all  of  the  kinetic  energy  of  the 
particle  is  dissipated  to  get 

Tzh  +  m^c^] 

E.=—, - ^  ^  ’ -  (3i 

Z\T  +  j  +  IbOOm^c^ 

where  is  the  energy  radiated  and  T  is  the  kinetic  energy  of  the  incident 
charged  particle.  If  we  assume  that  the  energy  radiated  during  deceleration  of 
the  charged  particle  is  radiated  as  a  single  coherent  photon,  then  (35)  gives  the 
energy  of  the  bremsstrahlung  photon. 
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The  origin  of  the  bremsstrahlung  photon  is  assnmed  to  be  at  the  point  of 


origin  of  the  decelerating  charged  particle  which  is  the  location  of  the  previons 
gamma  interaction.  This  assnmes  that  the  range  of  the  electron  in  the  scattering 
medinm  is  negligible  compared  to  the  ranges  of  the  incident  and  departing 
photons.  For  example,  an  electron  with  1  MeV  of  kinetic  energy  at  an  altitnde  of 
51  km  (where  the  air  density  of  the  US  Standard  Atmosphere  is  «  8.62  x  10”^ 
g/cm^)  and  electron  the  range-energy  relationship  for  air  from  NIST  (2010),  has  a 
total  path  length  on  the  order  of  1  km.  The  projected  range  of  the  electron,  or 
the  range  along  the  initial  direction  of  travel,  is  significantly  shorter  dne  to  the 
erratic  natnre  of  electron  paths  (Krane,  1988,  p.  197).  A  1  MeV  gamma  ray,  at 
the  same  altitnde,  has  an  average  mean  free  path  on  the  order  of  100  km.  Thns 
we  assnme  that  the  range  of  the  electron  is  insignificant  compared  to  that  of  the 
photon. 

Decay  of  Excited  Nuclei  by  Gamma  Emission 

The  emission  of  gamma  ray  photons  as  a  result  of  gamma  decay  of  excited 
nuclear  states  was  also  implemented  in  the  estimator  code.  Gamma  decay  level 
lists  were  acquired  from  the  Nudat  2.5  database  (NNDC,  2010)  and  converted  to 
a  format  to  be  read  into  the  estimator  code.  Level  lists  were  converted  to  a  text 
file  with  five  columns:  level  number,  level  energy,  emitted  gamma  energy, 
emission  intensity,  and  the  number  of  the  level  to  which  the  nucleus  decays  with 
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the  emission  of  that  gamma.  There  are  as  many  lines  in  the  file  for  each  level  as 


there  are  possible  gamma  emissions  at  that  level.  The  formatted  list  for  the  first 


six  levels  of  is  shown  in  Figure  19.  Note  the  emission  intensity  values  in  the 


fourth  column. 


Level 

Level  E(kev) 

Gamma  E(kev) 

Gamma  l 

New  Level 

1 

5270.155 

5269.161 

100.0 

0 

2 

5298. 822 

5297.817 

100.0 

0 

3 

6323.78 

1024.92 

0.05 

2 

3 

6323.78 

1053. 58 

0.1 

1 

3 

6323. 78 

6322. 35 

100.0 

0 

4 

7155.05 

831.27 

0.  5 

3 

4 

7155.05 

1856.11 

4.0 

2 

4 

7155.05 

1884.77 

100.0 

1 

4 

7155.05 

7153. 22 

0.023 

0 

5 

7300.83 

977.02 

0.25 

3 

5 

7300.83 

2001.86 

0.2 

2 

5 

7300. 83 

2030. 53 

0.6 

1 

5 

7300. 83 

7298.92 

100.0 

0 

6 

7567.1 

1243.2 

0.6 

3 

6 

7567.1 

2268.1 

4.0 

2 

6 

7567.1 

2296.8 

100.0 

1 

6 

7567.1 

7565.0 

1.  3 

0 

Figure  19.  Formatted  List  of  Levels  1-6  for 


Emission  intensities  for  each  emitted  gamma  at  a  level  are  given  relative  to 


the  most  likely  emitted  gamma  at  that  level,  which  has  emission  intensity  100. 


The  relative  distribution  of  emission  intensities  for  a  level  can  easily  be  converted 


to  a  probability  of  emission  by 


P 


7,2 


I  . 

7,2 


22, 


J  =  1 


(36) 


where  is  the  probability  of  emitting  the  zth  gamma  in  the  list,  is  the 
emission  intensity  of  the  zth  gamma  in  the  list,  is  the  number  of  gamma 
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emission  possibilities  for  the  level,  and  is  the  sum  of  all  emission 

j=i 

intensities  in  the  list  for  that  level. 

Inelastic  Neutron  Scatter  Gammas 

When  a  neutron  inelastic  scatter  event  is  simulated  by  the  estimator  code, 
the  material  and  inelastic  level  are  known.  We  then  use  this  information  to 
simulate  the  emission  of  photons  due  to  gamma  decay  of  the  excited  state  at  the 
site  of  the  inelastic  scatter  interaction. 

Although  gamma  emissions  by  gamma  decay  have  associated  half-lives 
before  emission,  we  assume  that  these  half-lives  are  insignificant  compared  to  the 
timescale  of  the  simulation  and  any  gamma  emissions  can  be  treated  as  occurring 
immediately.  This  assumption  is  appropriate  because  half-lives  for  gamma  decay 
are  generally  on  the  order  of  femtoseconds,  but  would  not  be  appropriate  if  any 
of  the  materials  included  in  the  simulation  have  metastable  states  with  half-lives 
significant  to  the  timescale  of  the  simulation.  We  also  assume  that  any  motion  of 
the  excited  nucleus  in  the  course  of  the  gamma  decay  is  negligible,  for  reasons 
similar  to  ignoring  the  path  length  of  charged  particles,  so  the  gamma  emissions 
are  simulated  at  the  site  of  the  inelastic  scatter. 

Using  the  known  material,  inelastic  level,  and  the  location  and  time  of  the 
inelastic  scatter  interaction,  the  estimator  code  draws  from  the  listed  gamma 
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emissions  for  that  level  according  to  their  probabilities.  The  selected  gamma 
emission  is  then  stored  for  later  transport  computations.  If  the  selected  gamma 
emission  decays  to  a  lower  excited  level  that  also  decays  by  gamma  emission,  an 
appropriate  gamma  emission  is  drawn  and  stored  again.  The  process  continues 
until  a  gamma  emission  is  drawn  that  decays  the  nucleus  to  the  ground  state, 
level  zero.  Alternatively,  to  eliminate  variance  from  sampling  the  gamma  decay 
emission  chain,  all  possible  gamma  emissions  weighted  according  to  their  relative 
likelihood  can  be  stored  for  later  computation. 

Neutron  Absorption 

Following  the  absorption  of  a  neutron  by  a  target  nucleus,  that  nucleus  may 
be  left  in  an  excited  state  that  decays  by  gamma  emission.  For  simplicity,  we 
exclude  the  radioactive  decay  of  an  unstable  isotope  created  by  neutron 
absorption,  and  include  only  the  decay  of  a  stable  nucleus  from  higher  energy 
states  by  gamma  emission.  For  implementation  in  the  estimator  code,  the  energy 
level  of  the  excited  nucleus  after  absorbing  a  neutron  is  the  energy  level  nearest 
to  but  not  greater  than  the  energy  of  the  absorbed  neutron  minus  the  Q-value  for 
the  neutron  absorption  interaction.  We  assume  that  this  residual  energy  is 
imparted  to  the  atomic  nucleus  by  placing  it  in  an  excited  state.  Recoil  of  the 
nucleus  and  other  energy  absorbing  processes  are  not  considered.  The  recoil 
energy  of  the  nucleus  in  neutron  absorption  interactions  is  not  negligible. 
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However,  we  ignore  it  here  to  simplify  the  physical  model  implemented  in  the 
code.  This  limits  the  physical  relevance  of  the  model,  bnt  implementation  in  the 
code  provides  the  infrastrnctnre  for  a  more  refined  model  to  be  implemented  in 
the  fntnre  Once  the  energy  level  of  the  excited  nncleus  is  determined,  gamma 
decay  can  be  modeled  similarly  to  gamma  decay  for  inelastic  scattering 
interactions,  bnt  nsing  the  level  list  for  the  isotope  of  the  absorber  atom  pins  one 
neutron. 

Gamma  Ray  Source  Energy  Selection 

Ray  Next-Event  Estimator  v4.0  included  several  distributions  for  drawing 
energy  of  a  particle  at  the  source  including  monoenergetic,  uniform.  Maxwell,  and 
Watt  distributions.  For  the  energy  of  gamma  rays  from  fission,  we  add  the 
following  distribution. 

A  fit  to  the  distribution  of  gamma  ray  energy  from  the  fission  of  is 
given  by  (Verbeke,  Hagmann,  &  Wright,  2010,  p.  12): 

38.13(e  -  0.085)e^®^^^  0.085  <  E  <  0.3 
A (e)  =  J  26.8e“^-^^  0.3  <  E  <  1.0  (37) 

8.0e'^-^^  1.0  <  E  <  8.0 

where  E  is  the  energy  of  the  gamma  ray  in  MeV.  This  distribution  can  be 
converted  to  a  probability  density  function  (PDF)  by 
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(38) 


=  4.547769900940718 
=  3.196439374382671 
=  0.954161007278409  . 

The  cumulative  probability  density  function  (CDF)  is  given  by 

(«> 

where  E'  is  a  dummy  variable  of  integration.  Selecting  integration  constants  so 
that  -P(-E')  remains  continuous  on  0.085  <E<  8.0  MeV ,  the  CDF  is 

C^+[C^E  0.085  <E<0.3 

P(e)  =  0.3<E<1.0  (41) 

1.0<E<8.0 

where 
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=  1.926282766794659 
a  =  2.759569114648494 

5 

a  =  1.909059196740567 

6 

=  0.850727368388764 
(7,  =  1.389756249731596 

o 

Cg  =  1.000130748747967 
=  0.867419097525827  . 

Figure  20  shows  plots  of  the  PDF  and  CDF  for  the  energy  distribution  of  gamma 
rays  from  the  fission  of 


Figure  20.  PDF  and  CDF  for  the  Distribution  of  Gamma  Ray  Energies  from  Fission  of 


To  sample  E  ,  we  draw  a  random  number,  ^  ,  uniformly  distributed 
between  zero  and  one  and  evaluate  P~^  to  find  E  .  Because  -P(-E')  is  a 
piecewise  function,  evaluation  of  must  be  treated  separately  for  each 
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interval.  For  0  <  ^  <  ,  or  random  nnmbers  corresponding  to  the  piece  of 

the  CDF  where  0.085  <  E  <  0.3  ,  P~^  can  be  evalnated  by  root-solving.  For 
P(0.3^  <  C  <  P^l.oj  and  P^l.oj  <  ^  <  1,  or  random  nnmbers  corresponding  to 
the  second  and  third  pieces  of  the  CDF,  P~^  can  be  evaluated  directly  by 


C7„ln 


P(0.3)  <  ^  <  i^(l.o) 

p(i.o)<e<i 


(42) 


where 


=  -0.434782608695652 
(7^2  =  0.719550640763897 
C^3  =  0.850727368388764 
(7^3  =  -0.909090909090909 
(7,,  =  1.152845265745635 

(7,«  =  1.000130748747967  . 

16 

Figure  21  shows  the  cumulative  distribution  of  100,000  energies  drawn  by  the 
estimator  code  from  the  distribution  in  (39). 

The  100,000  gamma  ray  energies  drawn  by  the  code  were  sorted  into  50 
equally  spaced  bins,  the  f^st  was  applied,  and  we  accepted  the  null 
hypothesis.  The  results  of  the  f^st  are  shown  in  Table  4. 
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P(E) 


Figure  21.  Cumulative  Distribution  of  100,000  Gamma  Ray  Energies  Drawn  by  the  Estimator 
Code  for  Fission  of  Compared  to  the  Gamma  Energy  CDF. 


Table  4.  Results  of  Chi-Squared  Test  for  Fit  of  100,000  Gamma  Ray  Energies  Drawn  by  the 

Estimator  Code  for  Fission  of 


x' 

Critical  Values  (90%  confidence,  49  DoF) 

Accept  Null 

Hypothesis? 

Upper 

Lower 

42.72 

66.34 

33.93 

Yes 

Summary  of  Major  Changes  to  RNEE  4.0 
Source  Code  Reorganization 

Because  RNEE  v4.0  was  developed  with  the  goal  of  demonstrating  the  ray 
next-event  estimation  method  for  neutron  transport,  significant  organizational 
changes  were  necessary  throughout  the  majority  of  the  source  code  to  allow 
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integration  of  transport  calculations  for  a  second  type  of  particle.  These 
organizational  changes  grouped  the  source  code  into  modules  by  function 
conducive  to  the  addition  of  gamma  ray  transport  computations.  The 
reorganization  also  served  to  create  an  open  program  architecture  to  allow  for 
future  modifications  and  additions  to  the  code. 

Variable  Dimensional  Extension 

Ray  Next-Event  Estimator  v4.0  used  several  large  multidimensional  arrays 
to  store  results  of  the  transport  computations.  These  allocatable  arrays  store  the 
computed  particle  flux  and  tally  count  for  each  time-energy  bin  at  the  flux-point 
for  each  batch  of  particle  histories.  To  utilize  the  existing  data  structure  as 
much  as  possible,  these  variables  were  extended  by  one  dimension.  Results  for 
neutron  transport  computations  were  stored  in  the  first  position  of  the  new 
dimension  and  results  for  gamma  transport  computations  were  stored  in  the 
second  position  of  the  new  dimension. 

However,  simply  adding  an  additional  dimension  index  would  have 
constrained  the  results  arrays  to  identical  dimensions  in  time  and  energy.  Energy 
range,  energy  resolution,  time  resolution,  and  time  of  interest  for  the  time-energy 
grid  at  the  flux-point  may  be  different  for  each  particle  type  included  in  the 
computation.  These  differences  result  in  different  time-energy  grid  dimensions 
for  each  particle  type  in  the  computation.  To  add  the  new  variable  dimension 
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and  preserve  the  flexibility  of  different  resnlts  array  dimensions,  a  Fortran 
derived  type  was  created  to  mirror  the  mnltidimensional  allocatable  array  in  the 
original  code.  A  new  variable  was  then  created  as  a  two-element  array  of  the 
new  type.  This  resulting  two-element  array  contained  two 

independently-allocatable  multidimensional  arrays  to  store  transport  computation 
results  for  the  two  different  particle  types. 

Photon  Interaction  Physics  Data 

Photo-Atomic  Cross  Sections 

Cross  sections  for  photo-atomic  interactions  were  acquired  from  the  XCOM 
Photon  Cross  Sections  Database  (NIST,  2010).  Interactions  included  are 
coherent  scattering,  incoherent  scattering,  photoelectric  absorption,  and  pair 
production  with  energy  ranging  from  1  keV  to  100  GeV.  The  cross  sections  can 
be  interpolated  using  logarithmic  functions.  The  cross  sections  are  queried  from 
the  online  database  and  downloadable  in  several  formats,  including  the  text 
format  used  by  RNEE  v5.1.  A  sample  is  shown  in  Figure  22. 
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Photon 

Energy 

l.OOOE-03 

1. 500E-03 

2.000E-03 

3.000E-03 

4.000E-03 

5.000E-03 

6.000E-03 

8.000E-03 

l.OOOE-02 

1. 500E-02 

2.000E-02 

3.000E-02 

4.000E-02 

5.000E-02 

6.000E-02 

8.000E-02 

l.OOOE-01 


coherent 
scatter . 
1. 29E+00 
1.18E+00 
1.05E+00 
8.00E-01 
6.11E-01 
4.77E-01 
3. 84E-01 
2. 69E-01 
2.03E-01 
1.21E-01 
8.04E-02 
4.23E-02 
2. 58E-02 
1. 74E-02 
1. 25E-02 
7. 30E-03 
4. 77E-03 


incoher . 
scatter. 
l.lOE-02 
2. 23E-02 
3. 51E-02 
5. 98E-02 
8.02E-02 
9. 57E-02 
1.07E-01 
1. 23E-01 
1. 33E-01 
1.48E-01 
1. 57E-01 
1. 63E-01 
1. 64E-01 
1. 62E-01 
1. 59E-01 
1. 53E-01 
1.46E-01 


Photoel . 
Absorb. 

3. 31E+03 
1.08E+03 
4. 76E+02 
1.45E+02 
6.10E+01 
3.09E+01 
1. 76E+01 
7.17E+00 
3. 54E+00 
9. 67E-01 
3. 81E-01 
l.OlE-01 
3. 91E-02 
1. 87E-02 
1.02E-02 
3. 92E-03 


Nuci ear 

Pr.  Prd. 

0. OOE+00 

0. OOE+00 

O.OOE+00 

O.OOE+00 

O.OOE+00 

0. OOE+00 

O.OOE+00 

O.OOE+00 

O.OOE+00 

0. OOE+00 

O.OOE+00 

O.OOE+00 

O.OOE+00 

0. OOE+00 

0. OOE+00 

0. OOE+00 

O.OOE+00 


Electron 
Pr.  Prd. 
0. OOE+00 
0. OOE+00 
0. OOE+00 
O.OOE+00 
O.OOE+00 
O.OOE+00 
O.OOE+00 
0. OOE+00 
O.OOE+00 
0. OOE+00 
O.OOE+00 
O.OOE+00 
0. OOE+00 
O.OOE+00 
O.OOE+00 
O.OOE+00 
O.OOE+00 


Tot.  w/ 

coherent 

3. 31E+03 

1.08E+03 

4. 77E+02 

1.46E+02 

6.17E+01 

3.14E+01 

1. 81E+01 

7. 56E+00 

3.88E+00 

1.24E+00 

6.18E-01 

3.07E-01 

2.29E-01 

1.98E-01 

1. 82E-01 

1.64E-01 

1. 53E-01 


Tot .  wo/ 
coherent 
3. 31E+03 
1.08E+03 
4. 76E+02 
1.45E+02 
6.10E+01 
3.10E+01 
1. 77E+01 
7. 29E+00 
3. 68E+00 
1.12E+00 
5. 37E-01 
2. 64E-01 
2.03E-01 
1. 81E-01 
1. 69E-01 
1. 57E-01 
1.48E-01 


1. 87E-03 

Figure  22.  Sample  of  Output  from  XCOM  Photon  Cross  Sections  Database 


Level  Lists  for  Gamma  Decay 

Lists  of  levels  for  gamma  decay  are  available  from  the  Nudat  2.5  database 
(NNDC,  2010).  The  data  is  available  for  each  isotope  in  tabular  html  format  or 
as  a  complete  evaluated  nuclear  structure  data  file  (ENSDF).  Because  only  a  few 
materials  were  used  in  the  course  of  this  research,  the  tabular  form  of  the  data 
was  converted  to  the  level  lists  used  to  implement  production  of  secondary 
gamma  rays  from  gamma  decay.  If  the  estimator  code  is  modified  for  application 
to  a  problem  with  many  atmospheric  constituents  or  requiring  level  data  for 
many  isotopes,  it  may  be  advantageous  to  modify  the  code  to  read  directly  from 
an  ENSDF  file. 
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Particle  Stacking  in  Linked  Lists 


Ray  Next-Event  Estimator  v4.0  only  considered  neutron  interactions  with  a 
single  particle  on  the  exit  channel.  Several  of  the  newly  added  gamma 
interactions,  as  well  as  the  inclusion  of  gamma  emission  cascades  from  neutron 
interactions,  have  the  potential  to  produce  multiple  photons  in  a  single 
interaction.  To  accommodate  these  groups  of  simultaneously  generated  particles, 
a  set  of  variables  was  added  to  the  code  to  serve  as  storage  for  particles  awaiting 
transport  calculations.  As  a  particle  is  produced,  either  at  the  source  or  in  a 
following  interaction,  it  is  placed  in  a  set  of  storage  locations  called  the  particle 
stack.  The  estimator  code  then  retrieves  a  particle  from  the  stack  and  proceeds 
to  make  the  transport  calculations  for  that  particle.  In  the  course  of  the 
transport  calculations,  more  particles  produced  in  interactions  of  the  current 
particle  may  be  placed  on  the  particle  stack.  After  completing  the  transport 
calculations  for  a  particle,  the  estimator  code  then  checks  for  another  particle  on 
the  stack  and  continues  the  transport  calculation.  If  no  particles  remain  on  the 
stack,  the  computation  is  complete. 

Memory  allocation  for  storage  of  particle  stacks  poses  an  interesting 
challenge  for  implementation:  How  much  stack  space  is  required?  The 
conditions  of  the  problem  being  solved  and  the  parameters  given  to  the  estimator 
code  can  cause  the  required  amount  of  storage  to  vary  wildly.  To  avoid 
quantifying  the  required  amount  of  storage  and  limiting  the  flexibility  of  the 


61 


estimator  code,  a  linked  list  was  implemented.  A  linked  list  is  a  dynamic  data 


strnctnre  in  which  each  element  in  the  list  contains  pointers  to  memory  locations 
of  the  previons  and  following  elements  in  the  list  (Chapman,  2008,  pp.  714-720). 
Linked  list  implementations  have  two  main  challenges: 

1.  The  process  of  incrementally  allocating  memory  for  each  variable  on 
the  stack  is  significantly  slower  than  the  single  large  array  stack 
implementation. 

2.  The  memory  nse  of  the  stack  is  nnbonnded  and  conld  become  too 
large  dnring  a  program  rnn. 

To  minimize  the  additional  time  reqnired  for  incremental  memory 
allocations,  the  management  of  the  linked  list  was  modified  so  that  when  a 
particle  was  read  from  the  list,  its  memory  was  not  deallocated.  Instead,  another 
pointer  variable  was  nsed  to  track  the  location  of  the  next  particle  to  be  read 
from  the  stack.  This  allowed  the  memory  allocated  for  the  stack  to  grow  and 
never  shrink.  When  a  new  particle  was  generated  to  be  placed  on  the  stack,  it  is 
now  stored  in  a  memory  location  that  is  already  allocated  and  the  pointer  to  the 
next  particle  to  be  read  from  the  stack  is  npdated.  Only  when  the  fnll  extent  of 
the  allocated  stack  memory  is  occnpied  by  particles  awaiting  transport 
compntations  will  a  new  particle  have  new  memory  allocated  for  it.  This 
minimizes  the  time  spent  allocating  memory  for  the  particle  stack  for  each  code 
rnn. 
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To  control  the  amount  of  memory  used  by  a  particle  stack,  counters  were 


implemented  to  track  the  number  of  elements  currently  in  the  list  and  the  total 
number  of  elements  allocated  for  the  list.  These  integers  could  then  be 
interrogated  to  determine  the  extent  of  the  particle  stack  memory  use  and  to 
force  memory  management  or  graceful  program  shutdown  should  the  memory 
allocation  for  the  list  become  excessive. 

Batch  and  Particle  Sequencing 

To  accommodate  the  coupled  neutron-photon  transport  problem  and  to 
allow  for  future  application  to  the  coupled  photon-neutron  problem,  neutron  and 
photon  transport  computations  were  implemented  in  a  phased  manner.  At  the 
start  of  a  run,  the  estimator  code  initializes  time-energy  grids  for  both  the 
neutron  and  photon  particle  types.  The  estimator  code  then  runs  the  neutron 
transport  problem,  if  one  is  specified.  During  the  course  of  the  computation,  if  a 
coupled  neutron-photon  problem  is  specified,  any  photons  induced  in  neutron 
interactions  are  also  transported  and  their  contributions  to  the  gamma  ray 
time-energy  grid  are  computed.  Once  computations  are  complete  for  the  neutron 
transport  problem  and  its  associated  gamma  ray  transport,  the  results  are 
processed  and  written  to  file. 

At  this  point,  if  no  additional  computations  are  specified,  the  program  run 
is  complete.  However,  if  a  gamma  ray  transport  problem  is  specified,  the 
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estimator  code  reinitializes  the  time-energy  grids  and  other  compntation  variables 


and  then  performs  the  gamma  ray  transport  compntation  in  a  similar  fashion  to 
the  nentron  transport  compntation  that  preceded  it.  In  this  manner,  the 
estimator  code  may  be  applied  to  conpled  nentron-photon  or  photon-nentron 
transport  compntations. 
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IV.  Estimator  Results  and  Demonstration 


The  Ray  Next-Event  Estimator  (RNEE)  v5.1  code  was  demonstrated  for 
neutron  transport  computations  and  compared  to  results  from  identical  problem 
runs  computed  by  RNEE  v4.0.  Both  code  versions  were  demonstrated  for 
transport  of  neutrons  from  a  14.06  MeV  monoenergetic  neutron  source  and  a 
fission  Watt  distributed  neutron  source.  After  determining  that  the  version  5.1 
was  functioning  consistently  with  version  4.0,  RNEE  v5.1  was  demonstrated  for 
gamma  ray  transport  from  a  15  MeV  monoenergetic  gamma  ray  source  and  a 
fission  distributed  gamma  ray  source.  Further  demonstration  was  made  for 
coupled  neutron-photon  computations  for  a  fission  Watt  distributed  neutron 
source  and  a  14.06  MeV  monoenergetic  neutron  source.  Table  5  lists  the 
properties  and  purpose  for  each  simulation  that  will  be  described  in  this  section. 

The  ray  next-event  estimator  code  produces  large  amounts  of  results  data 
including  time-energy  grid  flux  estimates,  time-energy  grid  contribution  tally 
counts,  statistical  comparisons  between  particle  batches,  and  time-energy  grid 
flux  from  direct  flight  contributions.  It  is  the  hope  of  the  author  that  the  limited 
data  presented  here  is  sufficient  to  illustrate  the  functionality  and  utility  of  the 
estimator  method  and  code. 
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Table  5.  List  of  Simulations  Used  to  Exercise  the  Estimator  Code. 


Run 

Code 

Version 

Description 

Pm-pose 

Total 

Histories 

Batches 

la 

4.0 

Neutron  Transport: 

Monoenergetic  Source 

v5.1  function  check 

2x10® 

40 

lb 

5.1 

2x10® 

40 

2a 

4.0 

Neutron  Transport: 

Fission  Source 

2x10® 

40 

2b 

5.1 

2x10® 

40 

3 

5.1 

Photon  Transport: 

Monoenergetic  Source 

Photon  Transport 

Demo 

1x10® 

20 

4 

5.1 

Photon  Transport: 

Fission  Source 

1x10® 

20 

5 

5.1 

Coupled  Transport: 

Fission  Neutron 

Source 

Coupled  Neutron- 

Photon  Transport 

Demo 

1x10® 

10 

6 

5.1 

Coupled  Transport: 

Monoenergetic 

Neutron  Source 

1x10® 

10 

RNEE  v5.1  Function  Check 

The  organizational  and  compntational  changes  to  the  estimator  code  were 
verified  by  two  test  simnlations:  Nentron  transport  from  a  monoenergetic  sonrce, 
and  nentron  transport  from  a  Watt  distribnted  fission  sonrce.  Each  simnlation 
was  rnn  on  version  4.0  and  version  5.1  of  the  estimator  code.  The  objective  was 
to  observe  near  identical  compntational  results  from  both  versions  of  the  code. 
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Run  #1  -  Neutron  Transport,  Monoenergetic  Source 


The  first  test  problem  was  the  transport  of  neutrons  from  a  14.06  MeV 
monoenergetic  neutron  point  source  at  50  km  altitude  to  a  flux-point  at  35,810 
km  altitude.  The  purpose  of  this  simulation  was  to  compare  the  results 
computed  by  version  5.1  of  the  estimator  code  to  result  of  an  identical  neutron 
transport  problem  run  from  version  4.0.  Properties  of  the  neutron  source  and 
flux-point  time-energy  grid  are  listed  in  Table  6.  Examples  of  input  files  to  both 
versions  of  the  estimator  can  be  found  in  Appendix  B. 


Table  6.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Monoenergetic 
Neutron  Source  Neutron  Transport  Demonstration. 


Source 

Neutron  Flux-Point  Time-Energy  Grid 

Type 

Neutron,  Monoenergetic 

Energy  Range 

0.01  -  15  MeV 

Energy 

14.06  MeV 

Energy  Resolution 

0.01  MeV 

^  of  Histories 

2  X  10®  in  40  batches 

Time  of  Interest 

20  sec 

Altitude 

50  km 

Time  Resolution 

0.01  sec 

Flux-Point  Altitude 

35,810  km 

The  scattered  free-field  flux  at  the  flux-point  per  source  neutron 
( n„  .  •  cm  •  cm”^  •  s~^  •  n‘^  1  for  each  time-energy  bin  is  plotted  for  versions 

\  ilux-pomt  source  J  ^ 

4.0  and  5.1  of  the  estimator  code  in  Figure  23  and  Figure  24  respectively. 
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20 

Figure  23.  Ray-Estimated  Neutron  Flux  Per  Source  Neutron  at  the  Flux-Point  from  a 
Monoenergetic  Point  Source  in  Time  and  Energy  Computed  by  RNEE  v4.0. 


20 

Figure  24.  Ray-Estimated  Neutron  Flux  Per  Source  Neutron  at  the  Flux-Point  from  a 
Monoenergetic  Point  Source  in  Time  and  Energy  Computed  by  RNEE  v5.1. 

By  graphical  inspection  and  nnmerical  comparison,  versions  4.0  and  5.1  of 
the  estimator  code  were  in  agreement  with  a  mean  relative  difference  of  0.66%. 
Perfect  agreement  was  not  expected  becanse  changes  to  the  code  strnctnre 
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changed  the  order  in  which  random  nnmbers  were  drawn  from  the  random 
nnmber  stream. 

The  ray-estimated  and  point-estimated  time  integrated  flnx  spectra 
compnted  by  code  version  5.1  are  plotted  Fignre  25  for  neutron  energies  up  to  2 
MeV  to  better  illustrate  the  behavior  of  neutron  contributions  to  the  flux  in  an 
energy  range  where  the  cross  section  varies  significantly  with  energy.  Note  the 
energy  correlated  peaks  in  cross  section  and  valleys  in  neutron  flux  at  the 
flux-point. 


«*cm 

Fluence  ( - ) 


cm  •  ^source 


Figure  25.  Ray-  and  Point-  Estimated  Time  Integrated  Flux  Spectra  for  Neutron  Transport 
from  a  Monoenergetic  Point  Source,  0  to  2  MeV. 
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Run  9^2  -  Neutron  Transport,  Fission  Source 


The  next  test  problem  was  the  same  as  run  ^1  except  the  neutrons  were 
emitted  from  a  Watt  distributed  fission  neutron  point  source.  The  purpose 
of  this  simulation  was  again  to  compare  the  results  computed  by  version  5.1  of 
the  estimator  code  to  result  of  an  identical  neutron  transport  problem  run  from 
version  4.0.  Properties  of  the  neutron  source  and  flux-point  time-energy  grid  are 
listed  in  Table  7. 


Table  7.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Fission  Neutron 
Source  Neutron  Transport  Demonstration. 


Source 

Neutron  Flux-Point  Time-Energy  Grid 

Type 

Neutron,  Fission 

Energy  Range 

0.01  -  15  MeV 

Watt  Distribution 

Energy  Resolution 

0.01  MeV 

of  Histories 

2  X  10®  in  40  batches 

Time  of  Interest 

20  sec 

Altitude 

50  km 

Time  Resolution 

0.01  sec 

Flux-Point  Altitude 

35,810  km 

The  Watt  distributed  fission  distribution  for  the  neutron  source  was 
built  in  to  version  4.0  of  the  estimator  code.  The  Watt  distribution  for 
fission  is: 

x[e)  =  0.453e“^“^®®  sinh (^|2.29E^  (43) 

where  E  is  neutron  energy  in  MeV  (Knief,  2008,  p.  45).  The  Watt  distribution 
is  plotted  in  Figure  26. 
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Figure  26.  Watt  Distribution  of  Neutron  Energies  from  Fission  of 


The  scattered  free-field  flux  at  the  flux-point  per  source  neutron 
( n„  .  •  cm  •  cm“^  •  s~^  •  n'^  1  for  each  time-enere;y  bin  is  plotted  for  versions 

\  ilux-pomt  source  J  ^ 

4.0  and  5.1  of  the  estimator  code  in  Figure  27  and  Figure  28  respectively. 

By  graphical  inspection  and  numerical  comparison,  versions  4.0  and  5.1  of 
the  estimator  code  were  in  agreement  with  a  mean  relative  difference  of  0.68%. 
Perfect  agreement  was  not  expected  because  changes  to  the  code  structure 
changed  the  order  in  which  random  numbers  were  drawn  from  the  random 
number  stream. 
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Figure  27.  Ray-Estimated  Neutron  Flux  Per  Source  Neutron  at  the  Flux-Point  from  a 
Watt  Distributed  Point  Source  in  Time  and  Energy  Computed  by  RNEE  v4.0. 


Figure  28.  Ray-Estimated  Neutron  Flux  Per  Source  Neutron  at  the  Flux-Point  from  a 
Watt  Distributed  Point  Source  in  Time  and  Energy  Computed  by  RNEE  v5.1. 


Based  on  comparison  of  estimator  resnlts  from  rnns  1  and  2,  and  many 
other  simnlations  rnn  in  the  conrse  of  this  research  and  thronghont  the  code 
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modification  process,  the  changes  between  estimator  code  versions  4.0  and  5.1  did 


not  impact  the  computational  results  of  the  code  for  neutron  transport. 

RNEE  v5.1  Demonstration:  Gamma  Ray  Transport 

Many  simulations  were  run  to  demonstrate  the  newly  implemented  gamma 
ray  transport  computation  capability.  Gamma  ray  transport  from  a 
monoenergetic  source  and  gamma  ray  transport  from  a  fission  source  are 
discussed  here. 

Run  9^3  -  Gamma  Ray  Transport,  Monoenergetic  Source 

The  first  gamma  ray  test  problem  was  the  transport  of  gamma  rays  from  a 
15  MeV  monoenergetic  gamma  ray  point  source  at  50  km  altitude  to  a  flux-point 
at  35,810  km  altitude.  The  purpose  of  this  simulation  was  to  demonstrate  the 
application  of  RNEE  v5.1  to  the  transport  of  primary  gamma  rays  from  a 
monoenergetic  point  source.  Properties  of  the  gamma  ray  source  and  flux-point 
time-energy  grid  are  listed  in  Table  8. 

The  scattered  free-field  flux  at  the  flux-point  per  source  gamma  ray 
( 7„  .  •  cm  •  cm~^  •  s~^  •  1  for  each  time-energy  bin  is  plotted  for  the  point 

\  '  imx-point  '  source  J  ^  ^ 

estimator  and  the  ray  estimator  in  Figure  29  and  Figure  30  respectively.  The 
inverse  relative  standard  deviation,  1  /  ,  for  each  time-energy  bin  is  plotted  in 
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Figure  31  to  compare  the  relative  certainty  of  the  flux  estimate  from  the  point 


and  ray  estimators  (larger  values  indicate  higher  certainty). 


Table  8.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Monoenergetic 
Gamma  Ray  Source  Gamma  Ray  Transport  Demonstration. 


Source 

Gamma  Flux-Point  Time-Energy  Grid 

Type 

Photon,  Monoenergetic 

Energy  Range 

0.01  -  15  MeV 

Energy 

15  MeV 

Energy  Resolution 

0.01  MeV 

^  of  Histories 

1x10®  in  20  batches 

Time  of  Interest 

0.2  sec 

Altitude 

50  km 

Time  Resolution 

0.00001  sec 

Flux-Point  Altitude 

35,810  km 

The  energy  peak  at  expected  from  annihilation  of  positrons  produced 

in  pair  production  interactions  was  not  observed  in  the  ray-estimated  flux.  One 
possible  explanation  for  the  lack  of  this  expected  feature  is  that  the  numerical 
quadrature  implemented  to  evaluate  the  ray  next-event  estimator  integrand  may 
not  appropriately  evaluate  the  endpoints  of  each  segment  of  the  ray.  These 
endpoints  are  explicitly  evaluated  in  the  point  estimator  by  definition,  giving  the 
point  estimate  the  distinct  energy  peak  at  . 
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Figure  29.  Point-Estimated  Gamma  Ray  Flux  Per  Source  Gamma  Ray  at  the  Flux-Point 
from  a  Monoenergetic  Point  Source  in  Time  and  Energy. 


Figure  30.  Ray-Estimated  Gamma  Ray  Flux  Per  Source  Gamma  Ray  at  the  Flux-Point  from 
a  Monoenergetic  Point  Source  in  Time  and  Energy. 


The  ray  estimate  also  showed  a  large  saddle  in  the  energy  range  1  to  5  MeV 
where  a  lower  flnx  at  the  flnx-point  was  estimated  than  in  the  point  estimate. 
The  variance  in  the  ray-estimated  flnx  was  also  mnch  higher  than  in  the 
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point-estimated  flux  in  this  region  of  the  spectra.  The  appearance  of  this  saddle 


may  be  an  artifact  of  the  angular  density  implementation  or  weighting  of  the 
likelihood  of  each  mechanism  as  the  estimate  is  computed  for  the  ray.  Because 
the  appearance  of  the  saddle  feature  is  limited  to  the  ray  estimate,  it  is  possible 
that  a  coding  implementation  for  one  or  more  photon  interaction  computations 
was  adequate  for  the  point  estimator,  but  not  the  ray  estimator. 


Figure  31.  Inverse  Relative  Standard  Deviation  of  Point-  and  Ray-  Estimated  Gamma  Ray 
Flux  at  the  Flux-Point  from  a  Monoenergetic  Point  Source  in  Time  and  Energy. 


We  expected  to  observe  lower  variance  in  the  estimate  computed  by  the  ray 
estimator.  However,  initial  investigation  indicated  that  the  variance  of  the  point 
estimate  was  generally  lower  along  the  ridge  of  peak  flux  and  in  other  areas  of 
the  grid  where  the  number  of  histories  contributing  to  a  particular  time-energy 
bin  was  high.  The  ray  estimator  provided  estimates  with  lower  variance  in 
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time-energy  bins  to  the  left  and  right  of  the  peak  flux  ridge  in  time  where  the 
number  of  contributing  histories  was  lower.  The  ray  estimator  had  lower 
variance  at  the  high  energy  end  of  the  spectrum  and  both  the  point  and  ray 
estimators  appeared  to  have  comparable  variance  at  energies  below  1  MeV. 

It  was  also  illustrative  to  compare  the  features  of  point  and  ray  estimates  as 
fluence  (flux  integrated  over  the  time  of  the  simulation).  The  ray-estimated  and 
point-estimated  time  integrated  flux  spectra  are  plotted  in  Figure  32.  The  energy 
peak  at  is  clearly  visible  in  the  point  estimate  and  missing  from  the  ray 

estimate.  The  low  energy  end  of  the  flux  saddle  is  also  visible  in  the  ray 
estimate. 
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Figure  32.  Ray-  and  Point-  Estimated  Time  Integrated  Flux  Spectra  for  Gamma  Ray 
Transport  from  a  Monoenergetic  Point  Source. 
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Run  9^4  -  Gamma  Ray  Transport,  Fission  Source 


The  next  test  problem  was  the  same  as  run  ^3  except  the  gamma  rays  were 
emitted  from  a  fission  gamma  ray  point  source.  The  purpose  of  this 
simulation  was  to  demonstrate  the  application  of  RNEE  v5.1  to  the  transport  of 
primary  gamma  rays  from  a  fission  distributed  point  source.  Properties  of  the 
gamma  ray  source  and  flux-point  time-energy  grid  are  listed  in  Table  9. 

The  scattered  free-field  flux  at  the  flux-point  per  source  gamma  ray 
( 7„  .  •  cm  •  cm“^  •  s"^  •  7“^  )  for  each  time-energy  bin  is  plotted  for  the  point 

\  '  rlux-point  '  source  J  ^  ^ 

estimator  and  the  ray  estimator  in  Figure  33  and  Figure  34  respectively.  The 
inverse  relative  standard  deviation,  1  /  ,  for  each  time-energy  bin  is  plotted  in 

Figure  35  to  compare  the  relative  certainty  of  the  flux  estimate  from  the  point 
and  ray  estimators  (larger  values  indicate  higher  certainty). 


Table  9.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Fission  Gamma 
Ray  Source  Gamma  Ray  Transport  Demonstration. 


Source 

Gamma  Flux-Point  Time-Energy  Grid 

Type 

Photon,  Fission 

Energy  Range 

0.01  -  8  MeV 

Energy  Resolution 

0.01  MeV 

^  of  Histories 

1x10®  in  20  batches 

Time  of  Interest 

0.2  sec 

Altitude 

50  km 

Time  Resolution 

0.00001  sec 

Flux-Point  Altitude 

35,810  km 
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Figure  33.  Point-Estimated  Gamma  Ray  Flux  Per  Source  Gamma  Ray  at  the  Flux-Point 
from  a  Fission  Point  Source  in  Time  and  Energy. 


Figure  34.  Ray-Estimated  Gamma  Ray  Flux  Per  Source  Gamma  Ray  at  the  Flux-Point  from 
a  Fission  Point  Source  in  Time  and  Energy. 
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Figure  35.  Inverse  Relative  Standard  Deviation  of  Point-  and  Ray-  Estimated  Gamma  Ray 
Flux  at  the  Flux-Point  from  a  Fission  Point  Source  in  Time  and  Energy. 

As  in  the  monoenergetic  gamma  ray  source  example,  we  observed  the 
missing  energy  peak  from  positron  annihilation  and  the  same  flux  saddle  in 

the  ray  estimate.  Without  comparison  to  an  estimate  computed  by  an 
independent  tool,  it  was  not  possible  to  determine  whether  the  point  estimator 
was  functioning  correctly  for  gamma  ray  transport,  but  the  results  of  the  point 
estimate  were  consistent  with  expectations  based  on  the  physics  of  the  transport 
problem.  If  we  assumed  the  point  estimator  was  working  properly,  then  these 
issues  with  the  ray  estimate  indicate  that  there  may  have  been  a  computational 
problem  with  the  behavior  of  photon  interaction  mechanics  as  the  ray  next-event 
estimator  integrand  was  evaluated  by  numerical  quadrature. 
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RNEE  v5.1  Demonstration:  Coupled  Neutron-Photon  Transport 


Rim  9^5  -  Coupled  Transport,  Fission  Nentron  Source 

The  first  coupled  neutron-photon  test  problem  was  the  transport  of 
secondary  gamma  rays  induced  by  prompt  neutrons  from  a  Watt  distributed 
fission  neutron  point  source.  The  purpose  of  this  simulation  was  to  demonstrate 
the  application  of  RNEE  v5.1  to  coupled  neutron-photon  transport  computations. 
Properties  of  the  neutron  source  and  flux-point  time-energy  grid  are  listed  in 
Table  10. 


Table  10.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Fission  Neutron 


Source  Coupled  Neutron-Photon  Transport  Demonstration. 


Source 

Gamma  Flux-Point  Time-Energy  Grid 

Type 

Neutron,  Fission 

Energy  Range 

0.01  -  8  MeV 

Watt  Distribution 

Energy  Resolution 

0.01  MeV 

of  Histories 

1  X  10®  in  10  batches 

Time  of  Interest 

1  sec 

Altitude 

50  km 

Time  Resolution 

0.0001  sec 

Flux-Point  Altitude 

35,810  km 

The  free-field  secondary  gamma  ray  flux  at  the  flux-point  per  source 
neutron  ( .  •  cm  •  cm"^  •  s~^  •  n"^  1  for  each  time-energy  bin  is  plotted  for 

\  '  liux-pomt  source  J  ^ 

the  point  estimator  and  the  ray  estimator  in  Figure  36  and  Figure  37 
respectively.  The  inverse  relative  standard  deviation,  1  /  ,  for  each 

time-energy  bin  is  plotted  in  Figure  38  to  compare  the  relative  certainty  of  the 
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flux  estimate  from  the  point  and  ray  estimators  (larger  values  indicate  higher 
certainty). 
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Figure  36.  Point-Estimated  Secondary  Gamma  Ray  Flux  Per  Source  Neutron  at  the 
Flux-Point  from  a  Fission  Neutron  Point  Source  in  Time  and  Energy. 
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Figure  37.  Ray-Estimated  Secondary  Gamma  Ray  Flux  Per  Source  Neutron  at  the 
Flux-Point  from  a  Fission  Neutron  Point  Source  in  Time  and  Energy. 
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Figure  38.  Inverse  Relative  Standard  Deviation  of  Point-  and  Ray-  Estimated  Secondary 
Gamma  Ray  Flux  at  the  Flux-Point  from  a  Fission  Neutron  Point  Source  in  Time  and 

Energy. 


It  was  immediately  apparent  that  both  the  point  and  ray  estimates  had  not 
converged  on  a  solntion  for  larger  values  of  time  and  energy.  Because  of  the  low 
mean  energy  of  neutrons  from  the  Watt  distribution,  the  1  x  10®  neutron  histories 
induced  a  total  of  only  0.24  x  10®  secondary  gamma  rays  for  transport 
computations.  The  low  number  of  secondary  gamma  rays  indicated  that  the  code 
implementation  for  secondary  gamma  ray  production  was  inefficient.  Longer 
computational  times,  for  corresponding  larger  numbers  of  prompt  neutron 
histories,  were  required  for  this  transport  problem. 


Run  9^6  -  Coupled  Transport,  Monoenergetic  Neutron  Source 


The  next  coupled  neutron-photon  test  problem  was  the  transport  of 
secondary  gamma  rays  induced  by  prompt  neutrons  from  a  monoenergetic 
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neutron  point  source.  The  purpose  of  this  simulation  was  to  supplement  the 
demonstration  of  run  ^5  with  a  coupled  neutron-photon  computation  with  a 
higher  number  of  induced  secondary  gamma  rays.  Properties  of  the  neutron 
source  and  flux-point  time-energy  grid  are  listed  in  Table  11. 


Table  11.  Properties  of  the  Source  and  Flux-Point  Time-Energy  Grid  for  Monenergetic 
Neutron  Source  Coupled  Neutron-Photon  Transport  Demonstration. 


Source 

Gamma  Flux-Point  Time-Energy  Grid 

Type 

Neutron,  Monoenergetic 

Energy  Range 

0.01  -  15  MeV 

Energy 

14.06  MeV 

Energy  Resolution 

0.01  MeV 

^  of  Histories 

1  X  10®  in  10  batches 

Time  of  Interest 

1  sec 

Altitude 

50  km 

Time  Resolution 

0.0001  sec 

Flux-Point  Altitude 

35,810  km 

The  free-field  secondary  gamma  ray  flux  at  the  flux-point  per  source 
neutron  ( .  •  cm  •  cm"^  •  s~^  •  n”^  1  for  each  time-energy  bin  is  plotted  for 

\  '  imx-pomt  source  J  ^ 

the  point  estimator  and  the  ray  estimator  in  Figure  39  and  Figure  40 
respectively.  The  inverse  relative  standard  deviation,  1  /  ,  for  each 

time-energy  bin  is  plotted  in  Figure  41  to  compare  the  relative  certainty  of  the 
flux  estimate  from  the  point  and  ray  estimators  (larger  values  indicate  higher 
certainty). 
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Figure  39.  Point-Estimated  Secondary  Gamma  Ray  Flux  Per  Source  Neutron  at  the 
Flux-Point  from  a  Monoenergetic  Neutron  Point  Source  in  Time  and  Energy. 


Figure  40.  Ray-Estimated  Secondary  Gamma  Ray  Flux  Per  Source  Neutron  at  the 
Flux-Point  from  a  Monoenergetic  Neutron  Point  Source  in  Time  and  Energy. 
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Figure  41.  Inverse  Relative  Standard  Deviation  of  Point-  and  Ray-  Estimated  Secondary 
Gamma  Ray  Flux  at  the  Flux-Point  from  a  Monoenergetic  Neutron  Point  Source  in  Time  and 

Energy. 

For  this  simulation,  the  1  x  10®  neutron  histories  produced  5.9  x  10® 
secondary  gamma  rays  for  transport  computations.  The  larger  number  of  gamma 
rays  produced  a  flux  estimate  with  much  higher  certainty  than  the  estimate  from 
run  ^5.  However,  the  variance  in  the  secondary  gamma  ray  flux  estimate  was 
still  higher  than  the  primary  gamma  ray  estimates  discussed  in  runs  3  and  4. 
Some  of  the  increase  in  variance  can  be  attributed  to  the  lower  number  of 
batches  used  to  compute  batch  statistics.  The  greater  variance  in  the  secondary 
gamma  estimates  is  also  a  product  of  the  fact  that  as  each  secondary  particle  is 
born  during  the  simulation,  it  inherits  the  weight  of  its  parent  at  the  time  of  its 
birth.  Since  the  parent  has  already  spent  some  simulation  time  transporting 
through  the  scattering  medium,  this  birth  weight  is  always  less  than  unity. 
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Thus,  despite  the  large  number  of  secondary  gamma  rays  produced,  each  of  these 


gamma  rays  carries  the  weight  of  less  than  one  whole  gamma  ray. 

Despite  the  higher  variance  in  the  estimated  flux,  many  energy  peaks  in  the 
secondary  gamma  ray  flux  estimate  were  apparent  in  the  point  estimate.  These 
peaks  correspond  to  gamma  emission  lines  from  the  decay  of  nuclei  excited  by 
neutron  interactions.  As  with  the  positron  annihilation  peak  in  the 

primary  gamma  ray  transport  demonstrations,  the  gamma  emission  lines  were 
not  visible  in  the  ray-estimated  secondary  gamma  ray  flux  estimate. 
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V.  Conclusions 


The  implementation  of  gamma  ray  interaction  physics  in  the  Ray 
Next- Event  Estimator  (RNEE)  v5.1  code  was  successful  with  limitations. 

Neutron  transport  computations  with  RNEE  v5.1  were  consistent  with  RNEE 
v4.0  despite  significant  organizational  changes  to  the  source  code.  Ray 
Next- Event  Estimator  v5.1  gamma  ray  transport  computations  using  the  built-in 
point  estimator  produced  flux  estimates  with  features  consistent  with 
expectations  based  on  the  physics  of  the  transport  problem.  Ray  Next-Event 
Estimator  v5.1  gamma  ray  transport  computations  using  the  ray  estimator  were 
not  consistent  with  expectations.  In  particular,  expected  energy  peaks  from 
positron  annihilation  and  gamma  decay  by  gamma  emission  were  not  observed  in 
the  ray  estimate.  The  ray-estimated  gamma  ray  flux  at  the  flux-point  due  to 
primary  gamma  rays,  runs  3  and  4,  also  included  the  low-flux  saddle  not  observed 
in  the  point  estimate.  Coupled  neutron-photon  transport  computations  were  also 
demonstrated,  with  the  observed  limitations  on  gamma  ray  flux  estimates  by  the 
ray  estimator,  and  a  need  for  reduced  computation  time  to  adequately  sample  the 
domain  for  secondary  gamma  ray  transport  computations  was  identified. 

Continued  development  of  this  estimation  method  is  important  because 
little  real-world  data  exists  with  which  to  validate  the  results  of  existing  or  new 
estimators  for  this  class  of  problems.  As  a  limited-utility  substitute  for  validation 
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against  real-world  data,  estimates  computed  using  different  methods  and/or  tools 
can  be  compared  with  one  another  for  investigation  of  the  utility  and  physical 
relevance  of  each  estimate.  With  further  refinement  of  the  physical  models  for 
neutron  and  photon  transport,  the  RNEE  code  will  become  a  formidable  tool  to 
be  used  in  conjunction  with  others  for  computational  transport  of  neutral 
particles  for  this  interesting  class  of  problems. 
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Appendix  A.  RNEE  5.1  Release  Notes 


These  release  notes  accompany  version  5.1.05  of  the  Ray  Next-Event  Estimator 
(RNEE)  Fortran  compnter  code.  Included  is  a  brief  description  of  the  program 
and  source  code,  known  issues,  and  recommended  areas  for  improvement.  This 
document  in  no  way  constitutes  a  user  manual  or  guide  to  the  use,  application, 
and/or  modification  of  the  RNEE. 


Ray  Next-Event  Estimator  v5.1.05 

Date:  06  February  2011 

Authors:  M.  Zickafoose,  W.  Dailey 

Advisor:  K.  Mathews 

Platform  of  Build:  Microsoft  Visual  Studio  2008  w/ 
Intel  Fortran  Compiler  12.0.2.154 
32-  and  64-  bit  architectmes 


The  source  code  is  comprised  of  31  modules  containing  subroutines,  functions, 
and  the  main  program  that  make  up  the  estimator  code.  The  ray-next-event 
estimator  has  been  developed  to  find  the  flux  at  a  specified  location  in  the  flux 
field,  which  we  have  called  the  flux-point.  In  the  code  implementation,  this  term 
is  awkward  as  a  variable  name.  Therefore,  the  term  detector  is  used  in  the  code 
and  file  formats,  and  consequently  in  these  appendices. 

RNEE  _  Driver 

•  main  program  module 

•  contains  version  details  and  history 

•  no  known  issues 

•  no  recommendations 

Monte  _  C  arlo  _  Module 

•  contains  top-level  routines  to  handle  running  of  simulation  batches  and 
managing  computational  data  flow 

•  no  known  issues 

•  Recommendation:  If  batch  or  particle  level  parallelism  is  implemented, 
most  of  the  main  parallel  statements  will  be  in  this  module. 


90 


Monte  _  Carlo  _  Neutron 

•  contains  top-level  routines  for  neutron  transport  computations 

•  module  variables  include  neutron  interaction  cross  sections 

•  no  known  issues 

•  Recommendation:  Should  additional  neutron  interaction  types  be  added, 
it  may  be  necessary  to  revisit  the  manner  in  which  neutron  cross  sections 
are  stored.  The  current  unified  energy  grid  is  a  simple,  logical  solution; 
but  it  is  memory  intensive  and  significant  computational  time  is  spent 
searching  the  cross  sections  tables. 

Monte  _  C  arlo  _  Gamma 

•  contains  top-level  routines  for  gamma  transport  computations 

•  module  variables  include  photo-atomic  cross  sections 

•  Known  Issue:  The  subroutine  Load_Cross_Sections  loads  the 
photo-atomic  cross  sections  into  a  unified  energy  grid  using  logarithmic 
interpolation  between  points.  Discontinuities  such  as  the  "k-edge"  are  not 
explicitly  handled:  The  loading  routine  does  not  look  for  special  cases. 

•  Recommendation:  Should  additional  photon  interactions  be  added, 
particularly  photo-nuclear  interactions,  the  loading  and  storage  of  cross 
sections  variables  will  need  to  be  reworked. 

Processing_Module 

•  contains  the  bulk  of  the  ray  next-event  estimator  computational  code 

•  computes  point  estimates,  ray  properties  and  geometry,  partitions  rays, 
and  computes  ray  estimates 

•  Parallel  estimation  is  implemented  in  the  subroutine 
P  art  it  ion_  and  _  Estimate 

•  no  known  issues 

•  Recommendation:  Review  the  implementation  of  Parallel  Estimation. 
Current  implementation  actually  degrades  performance  (longer 
computation  time)  if  enabled  for  simulations  with  low  resolution  detector 
grids. 

Sampling_Module 

•  Contains  most  of  the  routines  and  functions  that  utilize  random  numbers 
from  the  random  number  generator  stream 

•  Selection  of  particle  energies,  directions,  and  interaction 
material/mechanism  and  other  statistical  quantities  are  handled  here 
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•  no  known  issues 

•  no  recommendations 

Result  s  _  Module 

•  Contains  subroutines  to  condense,  compute,  and  write  output  data  to  file 

•  no  known  issues 

•  no  recommendations 

Transport  _  Functions  _  Shared 

•  Contains  subroutines  and  functions  shared  by  neutron  and  gamma 
transport  computations 

•  Includes  routines  with  code  specific  to  a  particular  particle  type,  but  each 
routine  is  used  to  transport  both  particle  types 

•  no  known  issues 

•  no  recommendations 

Transport  _  Functions  _  C  ommon 

•  Contains  subroutines  and  functions  common  to  neutron  and  gamma 
transport  computations 

•  No  particle  type  specific  code 

•  no  known  issues 

•  no  recommendations 

Transport  _  Functions  _  N  eutr  on 

•  Contains  subroutines  and  functions  specific  to  neutron  transport 

•  no  known  issues 

•  no  recommendations 

Transport  _  Functions  _  Gamma 

•  Contains  subroutines  and  functions  specific  to  gamma  ray  transport 

•  no  known  issues 

•  no  recommendations 

Root  _  Solver  _  Module 

•  Contains  functions  and  subroutines  for  root-solving 

•  no  know  issues 

•  no  recommendations 

Tallies  Module 
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•  Contains  rontines  to  initialize,  increment,  and  ontpnt  tallies  for  each 
detector  bin 

•  Modnle  variables  inclnde  tally  grids  for  the  point  estimator  and  the  ray 
estimator 

•  no  known  issues 

•  no  recommendations 


Setups 

•  Contains  setup  and  initialization  routines  for  the  estimator 

•  Reads  in  data  from  the  three  main  input  files:  ControlRNEE, 
Problem_Definition,  and  Method_Definition 

•  no  known  issues 

•  no  recommendations 

Statistics_Module 

•  Performs  statitistical  comparison  between  batches  and  writes  results  to  file 

•  no  known  issues 

•  no  recommendations 

Geometry  _  Functions 

•  Contains  low-level  functions  and  routines  for  computing  geometric 
quantities  for  Cartesian  and  spherical  geometries 

•  no  known  issues 

•  no  recommendations 

TE_  Spectra 

•  Contains  routines  to  initialize,  increment,  maintain  the  detector 
contribution  time-energy  spectra  grids 

•  Module  variables  include  batch-time-energy  spectra  for  each  detector 

•  no  known  issues 

•  Recommendation:  Implement  storage  and  maintenance  of  time-energy 
spectra  as  sparse  matrices.  This  will  reduce  the  memory  requirements  for 
the  program  and  improve  performance.  Also  look  into  the  need  to 
maintain  batch  results  in  memory.  Each  batch  could  be  written  to  file 
and  then  read  back  in  for  statistics  and  post-processing. 

Daughter  _  Particles  _  Module 

•  Contains  subroutines  to  sample  and  compute  secondary  particles  and  then 
store  that  particle  data  for  later  transport  computations 
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•  Known  Issue:  The  current  implementation  of  excited  state  gammas  uses  a 
recursive  subroutine  when  Level_Gamma_Decay_Sampled  is  set  to  False. 
If  an  atmospheric  constituent  is  added  with  many  excited  energy  levels, 
this  recursion  could  go  too  deep  and  cause  excessive  memory  use  or 
program  performance  degradation.  No  issues  were  observed  with  nitrogen 
or  oxygen. 

•  no  recommendations 

Table_  Search_  Module 

•  Contains  subroutines  and  functions  to  perform  searches  of  detector  grids 
and  cross  sections  variables 

•  no  known  issues 

•  Recommendation:  Although  an  efficient  binary  search  routine  forms  the 
backbone  of  table  searching  for  the  estimator,  the  majority  of 
computational  time  is  spent  searching  for  detector  time-energy  bins.  Any 
performance  enhancements  (maybe  parallel  searching?)  to  the  search 
routines  decrease  estimator  computational  time. 

Distribution_  Functions 

•  Contains  functions  to  compute  and/or  sample  from  various  energy 
distributions  for  initial  particle  energies 

•  no  known  issues 

•  Recommendation:  Renaming  of  distribution  functions  to  distinguish 
between  general  distributions,  neutron  specific  distributions,  and  gamma 
specific  distributions.  Also  refine  the  input  strings  in 
Problem_Definition.txt  where  the  source  distribution  is  specified. 

Error  _  Log  _  Module 

•  Contains  routines  to  write  error  conditions  to  file 

•  no  known  issues 

•  no  recommendations 

Quadratures  _  Module 

•  Contains  numerical  quadrature  routines,  general  purpose  and  customized 
for  RNEE  functions 

•  no  known  issues 

•  no  recommendations 

Particle  Stacks 
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•  Contains  rontines  to  initialize,  npdate,  and  track  memory  for  storage  of 
particle  awaiting  transport 

•  Modnle  variables  inclnde  particle  heaps  for  photons  and  nentrons 

•  no  known  issnes 

•  no  recommendations 

Interpolation_  Module 

•  Contains  fnnctions  to  linearly  and  logarithmically  interpolate  between 
data  points 

•  no  know  issues 

•  no  recommendations 

Phy  sical_  Constants 

•  Contains  parameters  for  physical  constant  values  and  program  run 
constants 

•  no  known  issues 

•  no  recommendations 

Prob_and_Meth_Vars 

•  Contains  variables  associated  with  the  problem  definition  including 
sources,  detectors,  and  geometry 

•  Contains  variables  associated  with  the  method  definition  including 
estimator  properties,  estimator  options,  and  particle/batch  specifications 

•  no  know  issues 

•  no  recommendations 

Types 

•  Contains  most  of  the  user  derived  types  used  in  the  program 

•  no  known  issues 

•  Recommendation:  A  critical  review  of  the  content  of  each  type  to 
eliminate  duplicated  and  unnecessary  data  structures.  This  could  help  to 
decrease  memory  requirements  and  improve  estimator  performance. 

Random_Numbers 

•  Contains  routines  and  functions  to  initialize,  maintain,  and  draw  from  the 
random  number  stream 

•  no  known  issues 
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•  Recommendations:  If  parallel  particles  or  batches  are  implemented,  these 
routine  will  need  to  be  modified  to  provide  independent  random  number 
streams  to  each  of  the  computing  threads. 

Kinds 

•  Contains  parameters  for  variable  types 

•  no  known  issues 

•  no  recommendations 
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Appendix  B.  Sample  Input  Files 


{This  page  intentionally  left  blank.) 
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RNEE  4.0  Control  File:  "Control_RNEE . txt" 
SRandomN umber Seed 
Seed_Input  =  T 


!  Can  be  True  of  False... Only  set  to  True  if  desire  correlated 
! sampling  between  RayNEE  and  PtNEE,  or  for  consistency  testing 
Seed  =  13  !  Use  a  prime  number  here,  if  the  above  value  is  T 

/ 

Slnputf lies 
num  of  materials  =  2 
num  of  cs  files  =  2 

num  of  cs  files  per  material  =  19,9,0 

num  of  geometries  =  2 

/ 

CO 
00 

SMaterialsList 

list  of  materials  =  'N-14', '0-16', '  ' 

A  value  of  materials  =  14 . DO, 16 . DO, 0 . DO 

material  fraction  =  0 . 788488D0, 0 . 211512D0, 0 . DO 

/ 

Slnputf lies 

!  Filename  for  the  number  density  integrals  --  Flat  Geometry 
! filename  =  'SlantIntNDbylayerFlat.txt' 
filename  =  'SlantINDbyHeightFlat.txt' 


!  must  be  same  size  array  as  the  list  of  materials 
land  num  of  materials,  and  in  same  order 
!  must  be  same  size  and  be  normalized  so  the  sum=l 


number  of  atmospheric  constituents  to  include 
includes  total,  elastic  and  level-inelastic  scatters 
number  of  files  for  each  constituent;  here  only  one 
constituent  but  19  cross  section  files 
loads  both  Flat  Geometry  and  Spherical  Geometry  Atmosphere 


/ 


Slnputf iles 

!  Filename  for  the  number  density  integrals  --  Spherical  Geometry 
! filename  =  'SlantIntNDbylayerSphr.txt' 
filename  =  'SlantINDbyHeightSphr.txt' 

/ 

SlnputFlles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\CrossSections\pendN14.txt' 

/ 

SlnputFiles 

CO 

^  !  Filename  for  the  material  and  mechanism 

filename  =  '.\CrossSections\pend016.txt' 

/ 

RNEE  4.0  Problem  Definition  File:  "  Problem  Definition.txt" 


SUnitsInfo 

one  per  cm2  =  T 

!  .TRUE,  puts  the 

output  in  units 

of  [ 1 /cm2 ] , 

.FALSE. 

leaves  it  in  units  of 

! [l/km2] 

counts  =  F 

/ 

!  .TRUE,  puts  output  in 

counts,  .FALSE. 

puts  it  in 

units  of 

Intensity 

SSourcelnf o 


100 


z  =  50 . DO 
X  =  O.DO 
y  =  O.DO 


[km]  -  geometric 
[km]  -  geometric 
[km]  -  geometric 

source  distribution  =  'Watt'  !  can  be  'Uniform',  'Watt',  'Maxwell',  'Point',  'Fusion',  'Group', 

! ' Boosted ' 

energy  =  14.06D  !  only  used  if  source  distribution  =  'Point' 

endpoints  =  1 1 . 7D0 , 1 1 . 7 IDO  !  only  used  if  source  distribution  =  'Group' 

sampled  direct  contribution  =  F  !  Used  to  determine  if  the  direct  flight  contribution  is  sampled 

/ 

!  The  following  Detector  namelist  is  for  a  spherical  geometry 
!  therefore  it  doesn't  contain  cartesian  coordinates, 

!  but  rather  the  height  or  look  angles... 

SDetectorlnf o 

z  =  35810. DO  !  [km] 

look  angle  =  44.8D0  !  [degrees] 

E_min  =  0.02D0  !  [MeV] 

E_max  =  15.D0  !  [MeV] 

T  max  =  2 O.DO  !  [s] 

E  spacing  =  'Linear'  !  Can  be  'Log'  or  'Linear' 

E_res  =  O.OIDO  !  [MeV]  only  used  in  E  spacing  is  'Linear' 

E  divisions  =  100  !  Integer  number  of  divisions  per  decade  starting  with  the  E  min  decade  and 

! ending  at  the  E  max  decade 

E  decades  =  4  !  Integer  number  of  decades  associated  with  the  E  min  and  E  max  values 

T_res  =  O.OIDO  !  [s] 

/ 


!  The  following  Detector  namelist  is  for  a  flat  geometry 

!  therefore  it  doesn't  contain  look  angles  and  the  height  is  set  negative 
!  but  rather  the  cartesian  coordinates... 

SDetectorlnf o 

[km]  --  but  negative  implies  Flat  geometry 
[km] 

[km] 

[km] 

[MeV] 

[MeV] 

[s] 

[MeV] 

[s] 


z  =  -1 .DO 
x_pos  =  O.DO 
y_pos  =  -35810. DO 
z_pos  =  35810. DO 
E  min  =  0 . 02D0 
E  max  =  15. DO 
T  max  =  2 O.DO 
E_res  =  O.OIDO 
T  res  =  O.OIDO 


/ 


RNEE  4.0  Method  Definition  File:  "Method  Definition.txt" 


SMethodDe tails 
Batch_size  =  100000 
num  of  batches  =  10 
interim  output  =  1000001 
diagnostics_basic  =  F 
diagnostics  full  =  F 

estimator  choice  =  'Both'  !  can  be  ' PtNEE ' , ' RayNEE ' ,  or  'Both' 

diagnostics  screen  =  T 
full  matrix  output  =  F 

Direct  Only  =  F  !  *Note  that  the  direct  contribution  is  computed  every  time,  but  this  flag 

!is  to  provide  a  sampled  direct  contribution  for  comparison  to  the 
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! computed  direct  contribution 

Integrand  Check  =  T 

/ 


SVarianceReduction 
Leakage  Suppression  =  T 
Absorption  Suppression  =  T 
Mat  mech  sampled  =  F 

/ 


!  F  ->  compute  all  materials  and  mechanisms  for  every  ray,  T  ->  sample  the 
Imaterial  and  mechanism  for  every  ray 


SDiagnosticsFileNames 

Basic  Ray  file  =  ' . \LogFile\Basic  Ray  Log.txt' 
Basic  Point  file  =  ' . \LogFile\Basic  Point  Log.txt' 
Full_File  =  ' . \LogFile\Full  Log.txt' 

/ 


RNEE  5.1  Control  File:  "Control_RNEE . txt" 

SRandomNumberSeed 
Use_MKL_RNGs  =  T 

Seed_Input  =  T  !  Can  be  True  of  False... Only  set  to  True  if  desire  correlated 

! sampling  between  RayNEE  and  PtNEE,  or  for  consistency  testing 
Seed  =  13  !  Use  a  prime  number  here,  if  the  above  value  is  T 

/ 


Slnputf lies 

num  of  materials  =  2 

num  of  neutron  cs  files  =  2 

num  of  gamma  cs  files  =  2 

num  of  cs  files  per  material  = 


number  of  atmospheric  constituents  to  include 

!  includes  total,  elastic  and  level-inelastic  scatters 
!  includes  total,  elastic  and  level-inelastic  scatters 
!  number  of  files  for  each  constituent;  here  only  one 


19,9,0 
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num  of  geometries  = 
/ 


! constituent  but  19  cross  section  files 
!  loads  both  Flat  Geometry  and  Spherical  Geometry  Atmosphere 
!  f ile  names 


&Mate rials List 

list_of  materials  =  'N-14', '0-16', '  ' 

A  value  of  materials  =  14 . DO, 16 . DO, 0 . DO  !  must  be  same  size  array  as  the  list  of  materials 

land  num  of  materials,  and  in  same  order 

material  fraction  =  0 . 788488D0, 0 . 211512D0, 0 . DO  !  must  be  same  size  and  be  normalized  so  the  sum=l 

/ 

SAtmosphereFiles 

!  Filename  for  the  number  density  integrals  --  Flat  Geometry 
filename  =  '.\AtmosphereFiles\SlantINDbyHeightFlat.txt' 

/ 

SAtmosphereFiles 

!  Filename  for  the  number  density  integrals  --  Spherical  Geometry 
filename  =  '.\AtmosphereFiles\SlantINDbyHeightSphr.txt' 

/ 

SNeutronCSFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\CrossSectionFiles\Neutron\pendN14.txt' 

/ 

SNeutronCSFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\CrossSectionFiles\Neutron\pend016.txt' 

/ 

SGammaCS Files 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\CrossSectionFiles\Gamma\N-atomic.txt' 

/ 


SGammaCS Files 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\CrossSectionFiles\Gamma\0-atomic.txt' 

/ 

SLevelFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\LevelsFiles\levelsN14.txt' 

/ 

SLevelFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\LevelsFiles\levels016.txt' 

/ 

SLevelFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\LevelsFiles\levelsN15.txt' 

/ 

SLevelFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  '.\LevelsFiles\levels017.txt' 

/ 

SLevelFiles 

!  Filename  for  the  material  and  mechanism 
filename  =  ' . \LevelsFiles\n  abs  Q.txt' 

/ 

RNEE  5.1  Problem  Definition  File:  "  Problem  Definition.txt" 

SUnitsInfo 

one  per  cm2  =  T  !  .TRUE,  puts  the  output  in  units  of  [l/cm2],  .FALSE,  leaves  it  in  units  of 

! [l/km2] 

counts  =  F  !  .TRUE,  puts  output  in  counts,  .FALSE,  puts  it  in  units  of  Intensity 

/ 
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&Neut ronS our ce Info 

z  =  50.D0  !  [km]  -  geometric 

X  =  O.DO  !  [km]  -  geometric 

y  =  O.DO  !  [km]  -  geometric 

source  distribution  =  'Point'  !  can  be  'Uniform',  'Watt',  'Maxwell',  'Point',  'Fusion',  'Group', 

! ' Boosted ' 

energy  =  14.06D  !  only  used  if  source_distribution  =  'Point' 

endpoints  =  1 1 . 7D0 , 1 1 . 7 IDO  !  only  used  if  source_distribution  =  'Group' 

/ 

SGammaS our ce Info 

!If  x,y,z  are  not  included  here,  they  will  be  the  same  as  the  previously  listed  source 
z  =  50.D0  !  [km]  -  geometric 

X  =  O.DO  !  [km]  -  geometric 

y  =  O.DO  !  [km]  -  geometric 

source  distribution  =  'Point'  !  can  be  'Uniform',  'Watt',  'Maxwell',  'Point',  'Fusion',  'Group', 

! ' Boosted ' , ' Fission! 3 5 ' 

energy  =  lO.ODO  !  only  used  if  source  distribution  =  'Point' 

endpoints  =  10 . ODO, 0 . OOIDO  !  only  used  if  source  distribution  =  'Group' 

/ 

!  The  following  Detector  namelist  is  for  a  spherical  geometry 
!  therefore  it  doesn't  contain  cartesian  coordinates, 

!  but  rather  the  height  or  look  angles... 

SNeutronDe tec tor Info 
z  =  35810. DO  !  [km] 

look  angle  =  44.8D0  !  [degrees] 

E  min  =  0.015D0  !  [MeV] 

E  max  =  15.D0  !  [MeV] 

T  max  =  2 O.DO  !  [s] 

E  spacing  =  'Linear'  !  Can  be  'Log'  or  'Linear' 

E  res  =  O.OIDO  !  [MeV]  only  used  in  E  spacing  is  'Linear' 

E  divisions  =  100  !  Integer  number  of  divisions  per  decade  starting  with  the  E  min  decade  and 

! ending  at  the  E  max  decade 

E  decades  =  4  !  Integer  number  of  decades  associated  with  the  E  min  and  E  max  values 

T_res  =  O.OIDO  !  [s] 

/ 
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!  The  following  Detector  namelist  is  for  a  flat  geometry 

!  therefore  it  doesn't  contain  look  angles  and  the  height  is  set  negative 
!  but  rather  the  cartesian  coordinates... 

SNeutronDe tec tor Info 

z  =  -l.DO  !  [km]  --  but  negative  implies  Flat  geometry 

X  pos  =  O.DO  !  [km] 

y_pos  =  -35810. DO  !  [km] 

z  pos  =  35810. DO  !  [km] 

E~min  =  0.02D0  !  [MeV] 

E_max  =  2 O.DO  !  [MeV] 

T  max  =  2 O.DO  !  [s] 

E~res  =  O.OIDO  !  [MeV] 

T_res  =  O.OIDO  !  [s] 

/ 

!  The  following  Detector  namelist  is  for  a  spherical  geometry 
!  therefore  it  doesn't  contain  cartesian  coordinates, 

!  but  rather  the  height  or  look  angles... 

SGammaDe tec tor Info 
z  =  35810. DO  !  [km] 

look  angle  =  44.8D0  !  [degrees] 

E_min  =  0.015D0  !  [MeV] 

E_max  =  15.D0  !  [MeV] 

T_max  =  l.ODO  !  [s] 

E  spacing  =  'Linear'  !  Can  be  'Log'  or  'Linear' 

E  res  =  O.OIDO  !  [MeV]  only  used  in  E  spacing  is  'Linear' 

E  divisions  =  100  !  Integer  number  of  divisions  per  decade  starting  with  the  E  min  decade  and 

! ending  at  the  E  max  decade 

E  decades  =  4  !  Integer  number  of  decades  associated  with  the  E  min  and  E  max  values 

T_res  =  O.OOOIDO  !  [s] 

/ 

!  The  following  Detector  namelist  is  for  a  flat  geometry 

!  therefore  it  doesn't  contain  look  angles  and  the  height  is  set  negative 
!  but  rather  the  cartesian  coordinates... 

SGammaDetectorlnf o 

z  =  -l.DO  !  [km]  --  but  negative  implies  Flat  geometry 

X  pos  =  O.DO  !  [km] 
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y  pos  = 

-35810. DO 

!  [km] 

z  pos  = 

35810. DO 

!  [km] 

E  min  = 

0.015D0 

!  [MeV] 

E  max  = 

15. DO 

!  [MeV] 

T  max  = 

0.5D0 

!  [s] 

E  res  = 

O.OIDO 

!  [MeV] 

T  res  = 

/ 

O.OOOIDO 

!  [s] 

RNEE  5.1  Method  Definition  File:  "Method  Definition.txt" 


SMethodDe tails 
geometry  =  "Sphr" 
diagnostics_basic  =  F 
diagnostics  full  =  F 
Integrand  Check  =  T 
estimator  choice  =  'Both' 
sampled  direct  contribution 
diagnostics  screen  =  T 
full  matrix  output  =  F 
Direct  Only  =  F 


secondary_particles  =  T 

Abs  Induced  Particles  =  T 
Level  Gamma  Decay_Sampled  = 

parallel  estimation  =  T 

/ 


!"Sphr"  or  "Flat"  selects  the  geometry  to  use  for  computations 


!  can  be  ' PtNEE ' , ' RayNEE ' ,  or  'Both' 

=  F  !  Used  to  determine  if  the  direct  flight  contribution  is  sampled 


!  *Note  that  the  direct  contribution  is  computed  every  time,  but  this  flag 
!is  to  provide  a  sampled  direct  contribution  for  comparison  to  the 
! computed  direct  contribution 

!  F  ->  suppresses  the  production  of  secondary  particles  in  all  reactions 
! (pair  production  still  produces  2  photons) 

!  must  have  secondary  particles  set  to  true  for  this  to  be  used 
F  !  F  ->  start  appropriately  weighted  photons  for  all  level  gamma  decay 
Ipossible  chains 


SVarianceReduction 
Leakage  Suppression  =  T 
Absorption  Suppression  =  T 

Mat  mech  sampled  =  F  !  F  ->  compute  all  materials  and  mechanisms  for  every  ray,  T  ->  sample  the 

Imaterial  and  mechanism  for  every  ray 


/ 
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SNeutronBatches 
Batch_size  =  100000 
num  of  batches  =  10 
interim  output  =  100001 
/ 

SGammaBatches 
Batch_size  =  100000 
num  of  batches  =  10 
interim  output  =  100001 
/ 

SDiagnosticsFileNames 

Basic  Ray  Diag  file  name  =  ' . \LogFile\Basic  Ray  Log.txt' 
Basic  Pt  Diag  file  name  =  ' . \LogFile\Basic  Point  Log.txt' 
Full  Diag  file  name  =  ' . \LogFile\Full  Log.txt' 

/ 
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